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Section i 


1.0 SUMMARY 

The main objective of the NASA-sponsored Aerothermal Modeling 
Program, Phase I was to assess current aerothermal submodels used 
in the Garrett Turbine Engine Company (GTEC) analytical combustor 
models . 

i'i number of "benchmark” quality test cases were selected 
after an extensive literature survey. The selected test cases, 
both nonreacting and reacting flows, were broadly divided into the 
following categories: 

o Simple flows 

o Complex nonswirling flows 

o Swirling flows 

o Dilution jet mixing in confined crossflows. 

These test cases were used to assess the following submodels 
separately and jointly for various combustion processes: 

o k-e model of turbulence and algebraic stress model, with 
and without various corrections including low Reynolds 
number and Richardson number corrections 

o Scalar transport models 

o Multistep kinetic schemes 

o Turbulence/ chemistry interaction 

o Spray combustion . 
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The following general conclusions were derived from Phase I 

work. 

o An accurate numerical scheme should be developed to mini- 
mize numerical diffusion in the computations of recircu- 
lating flows 

o Benchmark quality data should be generated under well- 
defined environments for validating the various sub- 
models used in gas turbine combustion analysis. 

o Although current aerothermal models make reasonable pre- 
dictions, intensive model development and validation 
effort should continue for the following submodels; 

Algebraic stress model 

Algebraic scalar transport model 

Two-step and four-step schemes 

Probability density function approach for a two- 
step scheme 

Double-reaction zone model. 
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SECTION II 


2,0 INTRODUCTION 

The objectives of the NASA Aerotherraal Modeling Program are to 
assess the current state-of-the-art and identify the deficiencies 
in current aerothermal models for gas turbine combustors. The pro- 
gram involves the following tasks: 

Task 1.1 - Model Definition 

Task 1.2 - Data Base Generation 

Task 1.3 - Benchmark Test Case Definition 

Task 2.1 - Model Execution 

Task 2.2 - Model Assessment 

Task 2.3 - Program Plan for Model Improvement. 

Paragraph 2.1 gives a brief background of aerothermal model- 
ing followed by a description of the Garrett empirical /analytical 
combustor design approach in Paragraph 2.2. This design approach 
is based on the use of a number of interrelated multidimensional 
analytical models that contain appropriate submodels (modules) of 
turbulence, chemistry, spray combustion/evaporation, soot, and 
high pressure radiation. These modules are described in Section 
3.0. A description of the numerical schemes employed are provided 
in Section 4.0, and a survey of relevant literature is presented in 
Section 5,0. 

The model assessment results are presented in four different 
sections: 

o Section 6.0 - Results for simple flows, with and without 
combustion 

o Section 7.0 - Results for complex nonswirling flows 
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o Section 8.0 - Evaluation of the models for swirling flows 

o Section 9.0 - Three-dimensional (3-D) dilution jet- 

mixing validation results 

Section 10.0 presents the conclusions and the recommendations for 
model improvements. 

2.1 Aerothermal Modeling Background 

Substantial increases in gas turbine performance have been 
achieved in recent years due largely to the use of advanced tech- 
nologies in components and material, in addition to operation at 
higher cycle pressures and temperatures. To meet the trend toward 
higher pressure ratio gas turbines with increased turbine inlet 
temperatures, increased research and development efforts have been 
directed toward the combustion system. These efforts have contri- 
buted largely toward gaining a better understanding of the overall 
combustion process and have led to the development of an advanced 
combustor design methodology based on a combination of empirical 
and analytical techniques. The challenging demands placed upon the 
combustion system due to increased performance and life require- 
ments, as well as the need to reduce combustor design and develop- 
ment cost, have provided the primary motivation for using multi- 
dimensional combustion analysis procedures. The advanced combus- 
tion analysis forms the basis for the design and development pro- 
cedures of advanced technology combustors at Garrett Turbine Engine 
Company » 

To provide greater confidence in the design of high-perform- 
ance, durable combustors for advanced aircraft turbine engines, a 
thorough understanding and accurate characterization of the various 
physical phenomena involved is required. Over the years, Garrett 
has been actively involved in the assessment, validation, and 


updating of combustor aerothermal models in the areas of multidi- 
mensional flow effects, effects of turbulence scale and intensity, 
combustion kinetics, fuel spray and flow field interactions, soot 
formation, and high-pressure flame radiation characteristics. 
Garrett has continued to assess every submodel within each model 
against fundamental data from ideal element tests. Concur- 
rently, model accuracy has been indirectly assessed by comparing 

predictions with measurements on a number of production and ad- 

1—9 

vanced combustors. Through an integrated effort of assessing 
both the models and the submodels, it has been possible to con- 
tinually improve the accuracy and reliability of the empirical/ 
analytical design procedure described in the following paragraph. 
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2. 2 Garrett Empirical/Analytical Combustor Design Approach 


Past approaches to the design and development of gas turbine 
combustion systems have largely involved the application of funda- 
mental knowledge of turbulent reacting flows on an empirical basis, 
followed by component testing to achieve optimum performance objec- 
tives. A number of semiempir ical relationships have been developed 
through the years to provide guidelines for the initial design of a 
new combustion system and to predict attainable performance on the 
basis of experience curves. Such an approach has been quite suc- 
cessful in the design and development of combustor configurations 
that are derived from proven concepts. 

The development of an empirical data base for combustors is 
evolutionary. Its limitations, regarding the development of 
advanced combustion systems with requirements outside of experience 
bounds, became apparent to Garrett in the early 1970's. The inade- 
quacy of the empirical approach in solving combustion development 
problems relating to gaseous and particulate emissions; carbon for- 
mation; and, more recently, liner and nozzle structural durability 
for high-temperature-rise applications required complementing this 
approach with advanced analytical methods. 

Garrett has developed a number of analytical models that form 
the basis for the design and development of advanced technology 
combustors. The internal flow field of modern gas turbine combus- 
tors is a highly complex 3-D phenomenon involving regions of 
reverse-flow. In addition, the various combustor regions require 
varying degrees of field resolution to predict accurately the con- 
vective and radiative fluxes. A modular approach, therefore, has 
been developed at Garrett allowing use of different computer 
models, as depicted in Figure 2.2-1. 
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ORIGINAL PAGE IS 
OF POOR QUALITY 


ANNULUS FLOW MODEL, AFM 


COMBUSTOR PERFORMANCE MODEL, CPM 

SOLVES FOR THE PRESSURE LOSSES 
AND AIRFLOW DISTRIBUTION WITHIN 
THE ANNULUS EXTERNAL TO THE 
COMBUSTOR. PREDICTS THE REQUIRED 
ORIFICE PATTERN FOR THE DESIRED 
FLOW SPLITS AND THE BOUNDARY 
CONDITIONS FOR THE COMBUSTOR 
PERFORMANCE MODEL. 


SOLVES THE GOVERNING REACTING 
FLUID DYNAMIC AND CHEMICAL REACTION 
EQUATIONS FOR THE ENTIRE COMBUSTOR. 
PREDICTS THE COMBUSTOR-FLOW FIELD 
INCLUDING VELOCITIES. TEMPERATURE, AND 
LOCAL FUEL-AIR RATIO AS WELL AS SMOKE 
AND SPECIES CONCENTRATION. 


/ 



TRANSITION MIXING MODEL, TMM 


SOLVES THE GOVERNING REACTING 
FLUID DYNAMIC AND CHEMICAL 
REACTION EQUATIONS IN THE 
180° TRANSITION LINER BEND 
PREDICTS THE TEMPERATURE 
QUALITY AT THE FIRST-STAGE 
TURBINE INLET. 


NEAR WALL MODEL, NWM 


SOLVES THE GOVERNING FLUID DYNAMIC 
AND HEAT FLUX EQUATIONS ADJACENT TO 
THE LINER WALLS USING A HluH-RESOLUTION 
GRID. PREDICTS LINER HEAT TRANSFER 
RATES AND ATTENDANT WALL TEMPERATURES. 


Figure 2.2-1 Combustor Models and Region of Application. 
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An annulus flow model is used to calculate pressure losses and 
airflow distribution within the annulus external to the combustor 
liner. This model calculates boundary conditions, such as flow 
distribution around the liner, jet velocity, and efflux angles, 
which are required as inputs for the combustor internal flow 
models . 

Two-dimensional (2-D) and 3-D combustor performance models are 
used to predict internal profiles of dependent variables including 
velocity, species, and temperature by solving fully coupled trans- 
port equations for turbulent, recirculating, spray-combusting flow 
fields. Up to 20,000 finite-difference grid nodes are numerically 
solved in these programs to ensure a relatively "grid-independent" 
solution for the main flow field. However, for the region close to 
the film-cooled wall, a better field resolution is required to 
accurately predict the convective fluxes and the wall temperatures. 
This is done by using near-wall models. 

The reverse-flow annular combustors generally employ transi- 
tion liners where the main flow direction changes from axial to 
radial for radial-inflow turbines or a full 18C-degree bend for 
axial flow turbines. The flow field has only small pockets of 
reverse-flow regions. Computationally more efficient 2-D and/or 
3-D transition mixing models are used for calculating the mixing 
rate of the cold dilution jets in the transition liner. These 
models calculate the turbine stator inlet profiles of temperature, 
velocity, and turbulence intensity, which are needed for assessing 
turbine hardware life. The various analytical models are used in 
the overall combustor design to arrive at a final combustor design 
in ja timely and cost-effective manner. 

The empirical/analytical combustor design approach is shown in 
Figure 2.2-2. The engine requirements and design define the com- 
bustor inlet conditions and limiting envelope constraints. Using 
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Figure 2.2-2. Coi..bustor Design Methodology. 














existing empirical design relations, a baseline combustor is 
defined and includes the appropriate flow splits, number of fuel 
nozzles, and orifice locations. The annulus flow model is then 
used to determine the orifice sizes to obtain the desired overall 
pressure drop, and flow splits needed to define the boundary con- 
ditions for the combustor performance model. 

The combustor performance model is run at various power condi- 
tions to evaluate combustor internal flow characteristics. If the 
design requires changes, the baseline combustor is altered, the 
annulus flow model is rerun, and the combustor performance model is 
again used to evaluate the new design. 

The combustor liner wall convective and radiative fluxes and 
attendant temperature levels and gradients are calculated with the 
near-wall model. The hot-side boundary conditions are defined by 
the combustor performance model. The cold-side boundary conditions 
are defined by the annulus flow model. The combustor performance 
model is also used to define initial conditions for the transition 
mixing model, which is used to calculate the mixing in the transi- 
tion liner and the resulting burner exit temperature quality. 

The results from the combustor performance model, near-wall 
model, and transition mixing model are factored into the analyti- 
cally designed combustor. If the design is lacking, iterations are 
peformed using the various models to arrive at an acceptable final 
configuration. 

The configuration is then fabricated and tested. If the 
result is unacceptable, the test data is compared with the analy- 
tical predictions and the appropriate subcomponent is modified, 
reanalyzed, and retested to verify that the modifications corrected 
the problems. This procedure is repeated until all the combustor 
design goals are achieved. Experience shows that this design 
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SECTION III 


3.C Description of Analytical Models 

Detailed descriptions of the combustor analytical models are 
provided in Paragraph 3.1. Each analytical model contains several 
submodels, which are described separately in Paragraphs 3.2 through 
3.6. 

3.1 The Analytical Models 

The Garrett modular analytical approach uses the following 
four models for analyzing gas turbine combustor flow field. 

o Annulus Flow Model (AFM) 

o Combustor Performance Model (CPM) , 2-D and 3-D 

o Transition Mixing Model (TMM) , 2-D and 3-D 

o Near-Wall Model (NWM) 

These models use submodels of turbulence, kinetics, radiation, 
and spray combustion/evaporation and dispersion as summarized in 
Table 1. 

3.1.1 Annulus Flow Model 


The first task in analyzing any combustion system is to pre- 
dict the annulus flow external to the combustor. For this, the AFM 
is used. The combustor annulus is divided into a number of sec- 
tions with the section boundaries defined by orifice rows in the 
liner or points of significant area change. In each section, the 
AFM solves the one-dimensional (1-D) equations for axial and tan- 
gential velocity. Mass is extracted from the annulus flow at each 
orifice row. The extracted mass is governed by the liner orifice 


PRECEDING PAGE BLANK. NOT FILMED 
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TABLE 1 


COMPUTER MODELS AND PHYSICAL SUBMODELS (MODULES) 


MODEL 

DIMEN 

SIONS 

NUMERICAL 

TYPE 

1 I 

SUBMODELS (MODULES) 

INPUTS 

INFO 

TO: 

OBTAINS 

INFO 

F p OM: 

TURBU 

LENCE 

KINETICS 

RADIA- 

TION 

1 

SPRAY 

ANNULUS 

1 


m , 

m 


_ 

CPM 

_ 

FLOW 









(AFM) 









COMBUSTOR 

2/3 

ELLIPTIC; 

K-e 

2-STEP/ 

4-FLUX/ 

LAGhANGIAN 

TMM 
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PERFORMANCE 
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• FLUX 
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geometry and semi-analytical correlations for discharge coeffi- 
cient. Pressure loss, both frictional and dump types, and heat 
transfer are included in the calculations. Iteration on combustor 
liner pressure drop continues until the total orifice flow rate 
achieves the desired value. The AFM predicts the static pressure 
distribution around the combustor, the liner pressure drop, orifice 
flow splits, and injection angles and velocities. These values are 
required as boundary conditions for the internal combustor flow 
programs (CPM, TMM, and NWM) . 

3.1.2 Combustor Performance Model 


The CPM has two versions: 2-D and 3“D. 
2-D Combustor Performance Model 


If the internal flow field of the combustor is predominantly 

2-D plane flow or axisymmetic flow, a 2-D CPM is used to calculate 

combustor internal flow field. The 2-D CPM is a generalized 

finite-difference program that solves the conservative form of the 

the governing fluid dynamic and chemical reaction equations, using 

15 

the numerical scheme of Patankar-Spalding. The following vari- 
ables are solved: 

o Axial, radial, and tangential velocity 

o Turbulent kinetic energy and dissipation 

o Total fuel, unburned fuel, and other chemical species 
including carton monoxide 

o Pressure 

o Stagnation enthalpy 
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o 


Radiation 


o Liquid particle droplet size, velocity, and evaporation 
rate 

Cylindrical or rectangular coordinates are used along with the 
capability of specifying any arbitrary shape for the liner wall or 
any arbitrary internal object such as a fuel nozzle shroud. Using 
the flow rates and velocities, etc., from the AFM, cooling slots, 
primary and dilution orifices, swirlers, and liquid or gaseous fuel 
nozzles are all modeled simultaneously. This gives the overall 
combustor flow field and the species and temperature distributions. 
Bulk flow properties determined include recirculation zone size and 
shape, primary and dilution jet penetration, and combustion effi- 
ciency. 

3-D Combustor Performance Model 

In many situations, the combustor geometry is not 2-D. In 

these cases, 3-D CPM must be used. The 3-D CPM is based on the 

12 

USARTL 3-D Model and can be considered an extension of the 2-D 
CPM to three dimensions. Both models use the same numerical scheme 
and the coordinate system. Like the 2-D CPM, the 3-D CPM solves for 
similar variables and requires boundary condition input from the 
AFM. Arbitrary complex boundaries and nozzle shrouds can be simu- 
lated. The 3-D CPM can analyze such 3-D flow situations as single 
(or multiple) swirlers in an annular combustor, tangential fuel 
nozzles, and discrete primary and dilution jets. 

3.1.3 Transition Mixing Model 

The TMM has two versions: 2-D and 3-D. 
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2- D Transition Mixing Models 

Though the 2-D CPM and 3-D CPM can analyze arbitrary shapes, 

i; ... 

they are limited to cylindrical or Cartesian coordinates. They 
cannot economically calculate the flow in the transition liner used 
in reverse- flow annular combustion systems. For relatively long 
combustors where the flow entering the transition liner is predom- 
inantly 2-D, the 2-D transition mixing model (2-D TMM) is used. 
This model is based on the GENMIX program of Patankar and 
Spalding^ . Modifications have been added that allow the program 
to negotiate 180-degree bends with the source terms added to 
account for the induced radial pressure gradients. Since it is a 
parabolic numerical scheme, this model is limited to transition 
liners in which the radii of curvature are large. Otherwise the 
pressure effects would have to propagate upstream. As in the other 
Garrett models, the two-equation k-e turbulence model is used along 
with the 2-step reaction mechanism. For initial profiles, the 2-D 
TMM uses the exit profiles as predicted by either the 2-D CPM or 3-D 
CPM. It then generates the exit profiles from the transition liner 
to which the turbine stator is exposed. 

3- D Transition Mixing Models 

With current trends toward shorter turbo-propulsion combus- 
tors and more compact transition liners, a significant amount of 
dilution jet mixing and spreading takes place within the transition 
liner. Attendant 3-D flow characteristics result from this mixing 
and spreading. Moreover, due to tight-bend radii of the transition 
liner, upstream (elliptic) effects caused by streamline curvature 
cannot be ignored. 

32 , A 3-D elliptic transition mixing model has therefore been 
developed that includes radiation and kinetic effects on the 
transition liner. This program is similar to the 3-D CPM, but 
has been structured to afford more than 2000 L finite-difference 
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nodes, where L is the number of nodes along the predominant flow 
direction. Theoretically, there is no limit on L; however, due to 
computer time consideration, L is generally kept less than 50. The 
3-D TMM can be adoptsd to analyze turbopropulsion combustors with 
much more complex geometries that cannot be adequately discretized 
by a cylindrical or Cartesian coordinate system. 

3.1.4 Near-Wall Models 

To accurately predict hot-side convective and radiative fluxes 
to the liner wall, a 2-D parabolic film cooling analytical model 
was developed during the Army Combustor Design Criteria Program. 
Subsequently, an improved 2-D NWM has been developed to allow a 
more accuitte assessment of the effects of the following on liner 
cooling effectiveness: 

o Slot geometry 
o Primary/dilution jets 
o Flow in the lateral directions 

o Radiation from the bulk flow field and the opposite wall 

o Spray combustion adjacent to the wall. 

The 2-D NWM can be used interactively with the combustor per- 
formance models to more accurately predict near-wall flow field. 
The 2-D NWM uses the same modules as the CPM. 

To further improve near-wall calculations, 2-D elliptic and 

l 

3-D parabolic NWM have also been developed at Garrett. 
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3.2 Description of Turbulence and Scalar Transport 


The internal flow field in gas turbine combustors is highly 
turbulent and recirculating. Efficient design of combustion 
systems requires a detailed understanding of the physio-cheinical 
processes of such systems. A prerequisite for this understanding 
is an ability to analyze the nonreacting turbulent recirculating 
flows. 

The fluid dynamics of turbulent flows are governed by the 
time-dependent Navier-Stokes equations. Solutions of these equa- 
tions are extremely difficult and require prohibitively large com- 
putational time. Furthermore, subgrid models are required to 
describe the transport phenomena in addition to the Navier-Stokes 
equations. A common alternative is to use time-averaged Navi.er- 
Stokes equations. This system of equations contains unknown higher 
order correlations resulting in a greater number of unknowns than 
the number of available equations. Turbulence models of the higher 
order correlations based on phenomenological assumptions are needed 
to close this system of equations. The degree of success of a tur- 
bulence model depends on the nature and accuracy of the phenomen- 
ological assumptions. 

The simplest of the turbulence models are the mixing-length 
models. In these models, the characteristic length scales of tur- 
bulence are often prescribed to close the system of equations. 
These models have been successful in treating simple flows like 
boundary layers and pipe flows, but have been unsuccessful in anal- 
yzing recirculating flows. 

The next higher order turbulence models are the one-equation 
models. These models solve one differential equation for deter- 
mining the distribution of the turbulent kinetic energy or equiva- 
lent characteristic property of turbulence. The characteristic 
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length scale is defined in a manner similar to mixing-length 
models. Consequently, the one-equation models have also been 
unsuccessful in predicting turbulent recirculating flows. 

The two-equation models are more complex than the mixing 
length models. These models use two differential equations for 
computing characteristic velocity and length scales. Among the 
two-equation models, the k-e model has been the most successful so 
far. The k-e model is used in the Garrett combustor analytical 
models and is described in Paragraph 3.2.1. In regions adjacent to 
walls, the viscous effects play a prominent role. To provide an 
accurate prediction of the flow in these regions, a low Reynolds 
number version of the k-e model is used in the Garrett near-wall 
models. This model is described in Paragraph 3.2.2. 

Even though the k-e model has been the most successful in pre- 
dicting recirculating flows, the predictions for flows with stream- 
line curvatures have been only qualitatively correct. Flow fields 
involving streamline curvatures have been known to have increased 
turbulence diffusion rates due to enhanced turbulence production. 
This increased turbulence production is not adequately accounted 
for in the k-e model. One way to include this extra production of 
turbulence is to modify the constants in the k-e model in propor- 
tion to the Richardson number, which is a measure of the extra 
strain rate produced by the streamline curvature. These correc- 
tions are described in Paragraph 3.2.3. 

The Richardson number corrections are applicable for flows 
with moderate streamline curvature effects. For flows with strong 
curvature effects, a solution of the Reynolds stress equations is 
necessary. Solution of the complete Reynolds stress components 
results in increased computational time. An attractive alternative 
is the use of an algebraic Reynolds stress model. In this model. 
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the terms in the Reynolds stress transport equations are approxi- 
mated to yield algebraic expressions in terms of the turbulence 
kinetic energy and length scale. The degree of approximations 
employed would determine the accuracy of predictions. The alge- 
braic Reynolds stress model is described in Paragraph 3.2.4. 

In complex combustor flow fields, the approximations used in 
developing the algebraic Reynolds stress model are not valid. For 
such flow fields, the Reynolds stress components must be obtained 
from the solution of differential transport equations for the 
appropriate Reynolds stress component. A description of these 
equations used by Garrett are provided in Paragraph 3.2.4. 

Another important submodel for combustor internal flows is the 
scalar transport model. The accuracy of the combustor performance 
predictions depends upon the accuracy of predicting the transport 
of scalar properties such as concentration of reactants, etc. The 
most commonly used scalar transport model is the gradient diffusion 
model. The gradient diffusion model does not adequately account 
for counter-gradient transport, which has been known to exist in 
most combustor internal flow fields. An algebraic scalar transport 
model (ASTM) has been developed at Garrett, which can predict 
counter-gradient scalar transport. This model is described in 
Paragraph 3.2.2. 

3.2.1 Governing Equations for the k-e Model 

The time-averaged transport equations for turbulence kinetic 
energy (k) and its dissipation rate (e) can be written in the fol- 
lowing generalized variable form: 
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Here p, r and denote the fluid density, the local effec- 
tive exchange coefficient of variable <£, and sources/sinks. The 
source terms for the dependent variables are 

I 

o Turbulence kintetic energy, k =^(u +v +w ) 

S k “ G k ~ p€ (2) 


where 
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o Dissipation rate. 


S, = (C, G, -c,pt)c 


The effective viscosity is obtained from the relation 


f £ eff = M+Mt 


where M and are the molecular and turbulent viscosities, respec- 
tively. is related to k and e via 

Mt=Copk2/ € (5) 

The exchange coefficients are defined as 

r eff,<£ ~ ^eff^°eff ,<£ 


Recommended values for the constants in the above equations are 


°eff 



= 0.09 

= 1.44 

= 1.92 

= 0.9 


^eff € ca l cu *st e d from 


where the K is the von Karman constant taken to be 0.41. 


( 6 ) 


The near-wall region is given a special treatment in the pro- 
gram. Since the expression for r ef f is accurate for turbulent 
flows only, a means is provided for the inclusion of the correct 
shear stresses and other fluxes at the wall. Therefore, the nodes 
next to the wall are assigned the following values as per an empi- 
rical wall law: 


^"• S Vwall =4 
y+>11 - 5 


Y* = P k l/2 C D l/A 4 


P ^ 9 ‘ 0( ^- ,,< ^ ) 

where 8 is the normal distance of the wall from the first interior 
adjacent node and cr is the laminar Schmidt number. 

The kinetic energy of turbulence has small diffusion near the 
wall; hence, r -q £° r ^ set equal to zero. Instead of computing 
r w?i , for c, it is calculated for the near-wall node by assuming a 
linear variation of the length scale giving the following expres- 
sion: 


e = C n 3/4 k 3/2 /(«5) 


3.2.2 Near-Wall Low Reynolds Number Correction 


In the near-wall region, the wall function approach, described 
in the previous paragraph does not properly describe the behavior 
of turbulence kinetic energy and its dissipation rate. A sys- 
tematic Taylor Series expansion technique has been developed by 
17 

Chien, which correctly describes the turbulent shear stresses, 
kinetic energy, and its dissipation rate in the region near a wall. 
To maintain the consistency of the behavior of k and e near the 
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wall, additional terms had to be added to the k and € equations and 
the turbulent effective diffusion rates were modified. The cor- 
rected equations of Chien near the wall are as follows; 

Turbulence kinetic energy; k 


TZ <<* k > + 7 ~§? <prV k) 




r 


d 

dr 


r r 


eff,k 


dk 

dr 


+ S k ♦ D 


where, D is the extra source term, given by 


D = 

y 

r eff,k = (lx + f » )/(r eff,k 
ffj = 1,0 -exp (-0.01 15 y + ) 


is the source term described in equation (2) 


Dissipation Rate, € 


■£ (fOO ♦ f -f- T (rPVe) = 
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eff,« ax 


♦i-S. 

r dr 


r Ti ae 

r eff,€ dr 


+ S + E 


(10) 


( 11 ) 


(12) 

(13) 


(14) 
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where, E is the additional source term for dissipation, defined as 


E = -2 /x exp (-0.5 y+) -4 

y 


(15) 


and 

with 


= ( C |G k “C 2 f 2 P^ 

(16) 

= 1.0 - 0.22 exp (-R t /6) 2 

(17) 

R P k2 

R T " /X£ 



In the modified k and £ equations, it is possible to apply the 
boundary conditions at the wall, with k = o and e = o at y = o. 
This approach gives consistent results near the wall. 

3.2.3 Richardson Numbers Correction to k-e Model 

The standard k- € model presented in Paragraph 3.2.1 describes 
the turbulence characteristics at any point in the flow field by a 
single velocity and length scale. These scales are obtained from 
an assumed isotropic turbulence structure. This model is adequate 
for simple flow fields. When significant streamline curvatures are 
introduced into the flow field, such as strong recirculation zones 
or swirl, the k-e model does not adequately account for the en- 
hanced turbulence diffusion caused by the extra strain rates asso- 
ciated with streamline curvature. For analyzing such flow fields, 
the k-€ model should be modified. 

A measure of the extra strain rate due to the streamline curv- 
ature is given by the Richardson number, Ri. The extra strain rate 
imposed on the flow field would tend to increase both the velocity 
and the length scales of turbulence. One way to account for the 
changes in the characteristic scales is to modify the turbulence 
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constants depending upon the Richardson number. For 2-D recir- 
culating flows, the Richardson number can be defined as 



k < R V 


(18) 


where, V R and R are the resultant velocity and the radius of curva 
ture respectively. They are defined as 


'r s \ f [ 


LT + V' 


(19) 


I 
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uv( |v _|y ) + U 2 |V 

9y ax Sx 


v 2 

. a. r 


( 20 ) 


For swirling flows, the corresponding Richardson number is defined 
by 

< V -f> Tr < rV »> 

' Ri v = (21) 

* ] ) 2 

In the k-e model, the governing equation for k is an exact 
equation and no empirical modeling is involved in it. However, the 
c equation is a modeled equation containing two empirical con- 
stants. The adopted approach is to modify the constant C 2 in the 
equation by the following expression 
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( 22 ) 


C 2 = 1. 92 exp (-« v Ri v - « Ri ). 

v d V e c c 

Here and a Q are empirical constants, whose values are of the 

order of 0.2. 

3*2.4 Algebraic Reynolds Stress Model 

Turbulent flow fields occurring in combustors are generally 
nonisotropic in character. The turbulent diffusion rates in dif- 
ferent directions are different depending upon the strain-rate ten- 
sor. Descriptions of such flow fields necessitate knowledge of the 
complete Reynolds stress components. Solution of the complete 
Reynolds stress components is expensive and complex. An alterna- 
tive to this approach is the algebraic Reynolds stress model. 

The algebraic Reynolds stress model is obtained by approximat- 
ing some of the higher order terms in the Reynolds stress equation 
based upon phenomenological simplifications. The approximations 
result in algebraic expressions for the Reynolds stress components 
with added empirical constants. 

The exact transport equation for Reynolds stress u~u7 at high 
Reynolds numbers in an incompressible flow can be written in the 
form 
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Pjj = Production 



7Tjj = Pressure-Strain 


( 23 ) 


At high Reynolds numbers, the viscous dissipation c i j is essential- 
ly due to a small scale turbulent motions and hence tends to iso- 
tropize u^u^ . The pressure-strain term has been modeled by Rodi 18 
in the form 





( 24 ) 


^ represents the contributions to pressure strain arising 
from fluctuating velocities only, and i r. . 0 accounts for the inter- 
actions between fluctuating velocities and mean strain. These 
terms are modeled as follows: 



-C' 
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Cu-u. 


3 6 ij k:1 



where: 
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Rod! has proposed that the convection and diffusion 
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With this simplification, the above equation reduces to the 
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(36) 


The ASM provides a mechanism by which anisotropic turbulence 
structure can be pred cted without substantial increase in computa- 
tional effort. The empirical constants in the ASM have not been 
fully established. They have to be determined by comparing the 
predictions with the data base. 
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3.2.5 Reynolds Stress Transport Model 


Although the algebraic Reynolds stress model provides a .means 
of computing the anisotropic structure of turbulence, its accuracy 
in complex turbulent flows is not expected to be good. In the alge- 
braic Reynolds stress model, the closure of the system of equations 
for Reynolds stresses was achieved by using a scaling law through 
which the higher order correlations were expressed as functions of 
lower order correlations. In complex internal combustor flows, the 
validity of the scaling laws is questionable. In st^ch cases, the 
only recourse available is to use the Reynolds stress transport 
model, where the Reynolds stress components are determined by solv- 
ing modeled differential transport equations for each stress compo- 
nent. The closure of these transport equations was achieved by 
modeling the higher order correlation terms in a manner analogous 
to the methods used in the k-e model. 


The governing equations for the Reynolds stress components can 
be written in generalized form as follows: 


(prW + -§g <PW# 


r | * ^ 

d ( r d$_\ _d_ ( r -p , JI&.) 

dx * rr eff, (}> dx' " dr Kr 1 eff,4> d r ' 

I _d_ /p 1 d4> J _ c 
' r dd u eff ,<f> r dd'\~ 4> 


(37) 


Here P, r eff and S 0 represent the fluid density, local effective 
exhange coefficient, and the source term for the dependent vari- 
able, <f>. The source terms are 

q Axial turbulence normal stress, u 2 
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Radial Turbulence normal stress, 2 
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+ Pvw (,-«-?) - c, 

Turbulence shear stress, uw 


vw 


s- = -P(l-o) 


+ P/3 


— au 2 av 0 — 

vw ir +u Sr +uv iF 


v. 


vw -w + * 2 air - uv v 


3V, 


Prk "ar - c 'i P c uw + 


(42) 


Here, C^, a, /3, and y are empirical constants, whose values are 
defined in Equation 26. 


Solutions obtained from these equations >are used in the mo- 
mentum equations instead of using the gradient diffusion assump- 
tions. In the new-wall region, boundary conditions for the depen- 
dent Reynolds stress components are applied by assuming the convec- 
tion and diffusion terms to be small, in accordance with the wall 
function approach. The Reynolds stress component w is computed 
from the relation 
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3.2.6 Scalar Transport Model 


The turbulent transport of Scalar properties in a flow is 
quite different from the transport of momenta. The most common 
method of describing turbulent scalar transport is the gradient 
transport law through the use of Prandtl/Schniidt numbers. This 
apprach is adopted in the standard k-e model. In the gradient 
transport model, the turbulent transport parameters of interest are 
defined by the following. 


where 
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*eff, 0 9x 

(43) 

pv0* 

se 

= " r eff,0 Sr 

(44) 

2 k 

dT\ 2 J 0 T\ 2 1 

(45) 

«0 € 

r eff, 0 [\ 8x / + \ 8r / J 



r ef.f ,0 = ^'eff/ <T 0 


Here, °0 is the Prandtl/Schmidt number for the scalar, 0 and 
a Q is an empirical constant. Recommended values for these are, 
°0 =0,9 and a ^ - 0.8 

The gradient diffusion model does not predict any counter- 
gradient diffusive transport. However, in many flows, regions of 
counter-gradient diffusive transport have been known to exist. For 
such flows, the scalar transport terms must be obtained from appro- 
priate transport equations. 
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3.2.7 Algebraic Scalar Transport Model 
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The governing equations for scalar transport are coupled non- 
linear differential equations, which are quite time consuming to 
solve. However, by using scaling laws analogous to the method 
used in developing the ASM, it is possible to obtain algebraic ex- 
pressions for the scalar transport correlations. These expressions 
could still account for the counter-gradient scalar transport. 
The detailed steps used in deriving the algebraic scalar transport 
terms are described in this paragraph. 
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a 1 and a 2 are empirical constants to be determined. 

1 2 

The transport equation for scalar fluctuations, 0 , is 


_ ,2 
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where , 
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(54) 


and 

e - a* ±- O' 2 

e 0" e k (55) 

For minimizing computational effort, the following assumption 
is made in a form consistent with the model for velocity fluctua- 
tions. 
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__ ,2 ~~a 

div (Pv0 ) - Diff (P 9 ) « C 0 (P- f ) ( 56) 

Using Equations (53) and (56) 

PFg " = ^0 ^ (P-e) (57) 

With these simplifications, the correlations for scalar quantities 
reduce to the form: 



Thr^ assumptions used in the derivation of the algebraic stress 
model are applicable for the flows that are close to local equi- 
librium. However, this model does not neglect any of the terms in 
the transport equation, and only a scaling law has been employed. 

The algebraic relations shown above express the turbulent scalar 
transport as a function of both mean scalar gradients and mean 
velocity gradients and hence can predict counter-gradient scalar 
transport. 
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3.3 Gaseous Fuel Combustion Models 


Successful modeling of combustors depends upon a correct de- 
scription and coupling of the fluid mechanics, turbulence, heat 
transfer, and chemical processes involved. The rates of turbulent 
exchange of various species and the rates of chemical change need 
to be modeled. In turn, this modeling will determine the details 
of quantities such as energy input to the gas stream, flow 
patterns, temperatures, and species concentrations. Turbulence 
models have been developed to a reasonably satisfactory stage. The 
state of development of chemical models is not nearly so satisfac- 
tory and is discussed in Paragraph 3.3.1. 

The turbulence/chemistry interaction model currently used by 
Garrett is a modified version of the eddy-breakup model; work is in 
progress on the development of a perturbation analysis technique. 
These models are described in Paragraph 3.3.2. 

3.3.1 Hydrocarbon Reaction Mechanisms 

A successful modeling of combustion systems depends on an 
adequate description of the reaction mechanism. For hydrocarbon 
oxidation, a large number of species participating simultaneously 
in numerous elementary kinetic steps is required to specify the 
reaction mechanism. These differential equations are ''stiff” and 
require special time-consuming integration methods. For a complex 
3-D problem, the computing costs would be prohibitive. Besides the 
large number of species equations to be solved, the elementary 
steps and their rate constants are not well known except for the 
simplest of hydrocarbons (e.g., CH^) . To avoid this problem, the 
gas turbine combustion modeling effort has frequently been simpli- 
fied by using a global approach that reduces chemistry to the 
specification of an overall global oxidation scheme. This can pre- 
dict quantities of interest: fuel consumption and heat release 
rates. 
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One Step Scheme 


The simplest global mechanism is the one-step scheme: 

C x H y + (X + |)o 2 -*xco 2 + \ h 2 o 

The advantage of this mechanism is its simplicity; it involves 
the solution of the conservation equations for unburned fuel and 
the mixture fraction. The heat release and the concentrations of 
the other species are then obtained from linear functions of the 
amount of fuel consumed. This mechanism, however, fails to predict 
the important characteristics of hydrocarbon oxidation, i.e., the 
formation of intermediates and CO, which influence the process con- 
siderably. As a result, this mechanism is inadequate for obtaining 
quantitative predictions. 

Two-Step Scheme 

A slightly more complex scheme is the two-step mechanism: 
C X H Y + <2 + 5) <°2 + nN 2 )—XCO + \ H 2 0 + (| + 5) nN 2 

X CO + £ ( 0 2 + nN 2 ) — X C 0 2 + ^ nN 2 - 

This scheme involves the solution of one additional equation: 
that for the concentration of CO. Although the two-step scheme has 
been widely used by Garrett, it is deficient in that the formation 
of intermediates is ignored. The derivation of the pertinent equa- 
tions is given below. 

For the first reaction, ' 
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f 


» (mass of 0 2 )/(mass of fuel) 
r 2 = (mass of CO)/(mass of fuel)* 

r 3 = (mass of H 2 0)/(mass of fuel); (61) 

in the second reaction: 

r 4 = (mass of o 2 )/(mass of CO) , 

= (mass of C0 2 )/(mass of CO). (62) 

The values of these ratios can be calculated in a straight-forward 
manner : 

r l “ (| + f> «0 2 /M fu 
r 2 - * "co^fu 
r 3 = <V2)W H20 /W fu 

r 4 = ( 1/ 2 ) Wq 2/^CO 
r 5 “ W C02^ W C0 

Here the W's are the molecular weights of the chemical species. 

The mass fractions of all chemical species obey the general 
differential equation with S^> as defined in Table 2. Further, the 
diffusion coefficient T ^ can be taken to be the same for all 
species, especially when the flow is turbulent. The value of is 
then M t /or t , where cr t is the turbulent Prandtl (or Schmidt) number. 
The source terms for various species are related via the ratios 
defined in Equations (61) and (62). As a result, the mass frac- 
tions of the species can be added in certain proportions to yield 


zero source terms. This is shown in Table 2. Here S fu denotes the 
source for fuel due to the first reaction, while S CQ stands for the 
rate of production of CO in the second reaction. 

TABLE 2. SOURCE TERMS FOR CHEMICAL SPECIES. 


<*> 

S 

M fu 

S fu 

m co 

S C0 “ r 2 s 

m ox 

r l S fu + r 4 ; 

M C02 

r 5 S C0 

m h 2 o 

r 3 S fu 



m ox 

- < ri + r 2 r 4 )m fu - r 4 m co 

0 



m C0 2 

+ r 5 m C0 + r 2 r 5 m £u 

\ 

0 


ill 

m H20 

+ r 3 m fu 

0 

K, 


The last three entries in Table 2 show that, because their 
source term is zero, a single solution for them would suffice pro- 
vided their boundary conditions are the same. . This condition can 
be ensured by normalizing the <£'s with reference to their values in 
the air and fuel streams. Thus a single variable f with a zero 
source term and with values 0 and 1 in the air and fuel streams 
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respectively can be regarded as providing the solutions for $ , <f» 

A B 

and <f> c via the following relationships: 


f = 




oir air 


r A r A,a r ' o J b, air 

^A, fuel" Air f B, fuel“^B,air 


^C,fuer^C,air 


(64) 


Further, let: 


^f^fuel 


= 1 


and 


ln Wair " R 


^2^ air 1_R 


(65) 


Combining Equations (64), (65), and the definitions of $ A , <f> B , $ c , 
we have: 


“ox = R(1 " £) + r 4 

m C0 + ^ r l +r 2 r 4^ 

< m £u- f > 

(66) 

m C0 2 m r 2 r 5 (£ -" 1 fu ) 

- r 5 m C0 


(67) 

"^2° = r 3 



(68) 

Incidentally, f can be 

considered as 

the mass fraction of 

" total 


fuel" that would prevail if the fuel did not react at all. 

The reaction rates S fjj and S CQ are given by the following 
relations: 

S fu = - (The smaller of and S 2 ) / 
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where 


s x - Fj_ P 1 ' 5 m ®- 5 m ox exp(-Ej/RT), 

5 2 * C R,1 "“fu </k - < 69 > 

Here, C R ^ is the eddy breakup constant for first reaction. 
Recommended value for C R ^ is 3.0 

S CO = - ( The smaller of S 3 and S 4 ), 

where 

5 3 = F 2 P m CO m OX ex P (~ E 2 //RT ^ ' 


S 4 = C p# 2 p m C0 ( 7 °) 

C R 2 is the eddy breakup constant for the second reation, recom- 
mended value for C R 2 is 4.0 


The constants in the 
values: * 

above expressions 

are 

given 

F 1 = 3.3xl0 14 , 

E^R = 27000., 

\ 

C R,1 

= 3, 

F 2 = 6 . 0 x 10 s , 
all in S.I. Units. 

E 2 /R = 12500., 

C R, 2 

11 


To summarize, the quantities m co , and f are used as the 
dependent variables of the differential equations. The source 
terms for m fu and m co are calculated from Equation (69) and (70), 
while the source for f is zero. The values of m 


‘OX 1 


m C02 ' an<3 m H20 


are then obtained from Equations (66) , (67) , and (68) . Lastly, m N2 
is calculated from the fact that all mass fractions should add up 
to unity. 
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Four-Step Scheme 


The oxidation of hydrocarbon fuel can be described by the fol- 
lowing steps: 

(a) Transformation of the hydrocarbon fuel into intermediate 
hydrocarbons and hydrogen with little release of energy 

(b) Oxidation of intermediates to CO and H 2 

(c) Oxidation of CO to CO 2 

(d) Oxidation of H 2 to H^O. 

Steps (b) through (d) are exothermic and are responsible for the 
release of energy and associated temperature rise. A reaction 
scheme, which is designed to correctly model the oxidation process, 
must include a description of these steps.. 

The simplest mechanism that accounts for the essential fea- 
tures of the hydrocarbon oxidation is the following four-step 

19 

scheme proposed by Hartman, et al. 


C N H 2N + 2“*I C 2 H 4 + H 2 
C 2 H 4 + 0 2 ~ 2C0 + 2H 2 

CO + 1/2 0 2 -> C0 2 

H 2 + 1/2 0 2 - H 2 Q 
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which is valid only for aliphatic hydrocarbons (of the type 
C N H 2N+2^ * To acGomn, °3ate a general hydrocarbon c x Hy„ the first 
two steps have been modified: 

C X H y C X H y- 2 + H 2 
C X H y -2 + l°2- XCO+ ¥ H 2 
This scheme involves the solution of two transport equations f 

V 

for the concentrations of and in addition to the trans- j 

port operations for unburned fuel, CO, and "total fuel" as outlined I 

in the two-step scheme. f 

3.3.2 Turbulence/Chemistry Interaction j 

In this section, a review of turbulent combustion models is | 

provided. This is followed by a description of the models under 
investigation at Garrett. Finally, a summary of turbulence/ 
chemistry interaction modeling is provided. 

Review of Turbulent Combustion Models 


An adequate treatment of turbulence/chemis|:ry interactions is 
essential for a reliable combustion model. S*ince the kinetic equa- 
tions are nonlinear in temperature and concentrations, large errors 
can be caused by incorrect time-averaging of the various terms with 
attendant effects on heat release rates, time-averaged gas temper- 
atures, and convective and radiative fluxes to the liner walls. 

The Problem - It has long been realized that the practice of writ- 
ing chemical kinetic equations in terms of time-averaged local var- 
iables such as 

W, = 7. 7. A exp | -0/T| 


(71) 


is unsound in turbulent mixing flows with relatively fast kinetics. 
Here VT is the chemical reaction rate for species i of mole frac- 
tion Y ^ ; Yj is the mole fraction of another species; A is the pre- 
exponential factor in the Arrhenius expression for the chemical 
kinetic rate; Q is the activation temperature; and V i % the abso- 
lute temperature, the overbar indicating time averaging. Equation 
(71) neglects the correlations between fluctuations in the various 
quantities, e.g., Y* jY' j ' Y t Tt t and the contributions from these 
terms can change the computed reaction rate by an order of magni- 
tude or more. Attempts to compute the various correlations direct- 
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iy, as has been done by Donaldson and Borghi, have proved suc- 
cessful only in flows where the fluctuations are low and the heat 
release is not large. 

The Fast-Chemistry Approach - A more satisfactory approach in non- 
premixed combustion systems is based on the assumption that the 
chemistry is fast. The chemical reaction rates are then entirely 
mixing controlled and are a function of the turbulence rather than 
the chemical kinetics. Two versions of this approach are in cur- 
rent use. 

In the first version, equations for conserved scalars such as 
the element mass fractions or the mixture fraction are solved, 
instead of solving directly for the individual combustion product 
species. Molecular species concentrations and temperature are then 
determined from the computed moments of the conserved scalar, 
usually by assuming some probability density function for the con- 
served scalar. The fast chemistry assumption implies an instanta- 
neous relationship between molecular species or temperature and the 
conserved scalar. Chemical reaction rates can be found, if needed, 

by differentiation. Such an approach has !??en used by some inves- 

22 , 

tigators, e.g., Lockwood. The prohxerr is that the extension to 
include the treatment of any tvp'a- reaction mechanism, even a 
single-step one, entails complication*?, * 
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In the second version, molecular equations are solved and the 
chemical reaction rate is modeled directly in terms of turbulence 
quantities. The Spalding eddy breakup (EBU) model is the prime 
example of this approach, in which single or multi-step reaction 
mechanisms can be handled, when suitable modifications to the model 
are made. This has been done by Garrett, which has developed a mod- 
ified version of this model. 

Both of these approaches give qualitatively satisfactory 
results for the main species concentrations and temperature. The 
problem, however, is that chemical kinetics is no longer involved. 

The Inclusion of Chemical Kinetics - Although the majority of fuel 
oxidation in gas turbine combustion systems is essentially mixing 
controlled under most operating conditions, the chemical kinetic 
effects must be included to predict hydrocarbon emissions, flame 
stabilization or blowout, CO emissions, soot formation and burnout, 
and NO formation. The problems of satisfactorily including the 
chemical kinetics into the chemical reaction rate have proved to be 
formidable. 

As a first step towards inclusion of the kinetics, the EBU 
model has been modified at Garrett to compute tke reaction rate R 
from the minimum of the EBU rate and the well-stirred reactor 
global reaction rates. Garrett has used 2-step and, recently, 4- 
step kinetic schemes. The procedure is illustrated 'here with a 
single-step reaction scheme. 


t 

~ m ' n £ R EBU’ r ws^ 

(72) 

‘EBU 

= Cj^P</>e/k 

(73) 

WS 

= A p' ' S M ox M fu 0 ’ 5 exp (-E/RT) 

(74) 


= min C M QX /i 3 

(75) 
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where: 


R EBU * 


HwS " 



A 

P 


M 


OX 


M fu 

E/R 

T 

k 


€ 


i 


the eddy-breakup rate of chemical reaction; 
the well-stirred rate of chemical reaction; 
empirical constant; 

Arrhenius pre- exponential factor; 
density; 

mass fraction of oxidant 

mass fraction of fuel 

activation temperature; 

absolute temperature; 

turbulence kinetic energy 

dissipation rate of k 

mass of oxidant per unit mass of fuel. 


This model, which was used in the computer codes developed by 

12 

Garrett for the US Army, has been found suitable for qualitative 

correlations. A further extension of the model at Garrett solves a 

2 

transport equation for the fluctuation, g (=0 ) , of the fuel con- 
centration rather than obtaining 0 from Equation (75) . This addi- 
tional equation results in better agreement between the predictions 
and experimental data for typical combustor geometries. It is 
still not suitable for accurate quantitative predictions and for 
problems such as kinetic effects on temperature or satisfactory 
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estimates of free-radical concentrations. These quantities are 
required for accurate prediction of CO, since CO consumption often 
occurs via the reaction CO + OH = C0 2 + H. 

Achievement of a satisfactory approach to the modeling of the 

chemical reactions has led to several rather novel approaches to 

the problem. Many of these methods, such as Spalding's ESCIMO 
24 

model, Chorin's Vortex Dynamics, large eddy simulation tech- 
niques, and joint PDF approaches, 3 involve at least an order of 
magnitude increase in the size and complexity of the computation 
and as yet are not completely formulated. One approach, based on a 
perturbation analysis for reaction kinetics, does not involve such 
an increase in size and complexity, and has been formulated by 
Bilger. This method, as described in the next paragraph, has 
been adapted by Garrett. 

The Perturbation Analysis for Reaction Kinetics 

In this approach, a turbulent diffusion flame model has been 
27 

developed. It uses a double reaction zone model and perturbation 
analysis for finite rate kinetics for hydrocarbon combustion. A 
system of eight parabolic transport equations is solved. The system 
consists of the usual k-e model equation in Favre averaged form for 
continuity, momentum, mean mixture fraction, specific turbulent 
kinetic energy, and turbulent kinetic energy dissipation rate, with 
additional scalar transport equations for mixture mole number per- 
tubation, unburned fuel mass fraction pertubafion, and mixture 
fraction variance. The thermodynamic state (and composition) of 
the flow field is contained in an equilibrium model of the hydro- 
carbon-air mixture in terms of mean mixture fraction. The progress 
of the chemical reactions (and thereby the molecular kinetic rates) 
is contained in perturbations or constraints on the equilibrium in 
terms of mole number and unburned fuel mass fraction. The unburned 
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fuel mass fraction and the intermediate are specified as functions 
of mixture mass fraction, £. For the fuel. 


. where 


and 


where 


Y 1 3 “ Y ° 1 3 + y 


(76) 


Y °I3 = 0 


0<£<£ig 


(f-U 

VO _ JS_ 

13 '^7 


^ig < £ < I 


y is the pertubation in fuel mass fraction 

Y° 13 is the "fast chemistry" (i.e., zero pertubation) fuel 
mass fraction 


Y 13 is the fuel mass fraction 
P . 

ig is the mixture fraction where a "reaction sheet" 
consumption or pyrolysis of fuel occurs under fast 
chemistry conditions similar to the classical 
Burke-Schumann formulation. here is taken as 

0.073. 


For the intermediate hydrocarbon, its mass fraction is given by. 



****** 


(77) 


Y 


12 = 


er 




*s<***lg 


(78) 



er ^rg-^ 

rrrgrhy 


<i -£) 


V 


£ ^1 


(79) 


where 

r is the mass fraction of fuel in the inlet fuel 
stream 

Y ^2 is the intermediate hydrocarbon mass (fraction 

e is the fraction of fuel by mass that is converted to 

intermediate at for fast chemistry conditions, 

e is taken as 0.2 


and 


£ s is the stoichiometric mixture fraction 


Hence, the double reaction zone at £ and £. . 

s ig 

Thirteen species are considered in the reactions. H, H 2 , H 2 0, 
0, OH, 0 2 , H0 2 , N«i , Ar, CO, and C0 2 are calculated from partial or 
constrained equilibrium, and the fuel and intermediate are speci- 
fied as in Equations 76-79. The pertubation in mole number space 
is a result of the rate of three-body recombination reactions. 







H+H+M 


H 2 + M 

R1 

H+OH+M 

— 

H 2 0 + M 

R2 

h+o 2 +m 

— ► 

H0 2 + M 

R3 

H+O+M 

r 

OH + M 

R4 


The progress of these reactions toward equilibrium is measured by 
mole number, N 


N = 


II Y : 


S 

1=1 


W. 


Mole 

“W 


(80) 


where 


and 


is the mass fraction of species i 


is the molecular weight of species i 


The pertubation in mole number is defined as 


Where N° is the number of moles at full equilibrium. Then, for vari- 
ous pertubations in mole number space for a given mixture mass 

fraction, the time rate of change of N can be calculated from reac- 
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tions R1-R4 using here the kinetic data of Jensen and Jones. 

Perturbation Equations 

From the species balance transport equation. 


aY ; 


at 


+ pU k 




= 0 


(82) 


where the molecular diffusivity is assumed to be the same for all 
species, and using Equations 80 and 81 gives. 


dn 

a t 


nl i dn 

pU k di 





(83) 


i 

w n is the source term for mole number 


A similar equation can be written for the fuel mass fraction, per- 
turbation. 




ay a 
ax k - ax k 





+ W. 


(84) 


Wy is the source term for mass fraction 
where the dependence of y on N° is neglected. 
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-rtsss* 

Then using density weighting (Favre averaging) , the equations for 
turbulent flow become. 



Using gradient modeling of the turbulent fluxes and k-e modeling of 
the scalar dissipation X 

where 



(87) 


the Favre averaged scalar dissipation rate is 

x -% H 


( 88 ) 


where C is a model constant with value of 1.79, Then, 


P u. 


an a 

k Sx k ~ ax k 


^eff 
n ax k 


an 

cr ax, 


1/2 P (c/k) . N>s 


(89) 


+ w 


p “ k s *k 


ax, 


^eff 


cr. 


SX 


= 1/2 P C 


92 


C" 2 (<A) Ay P(f g ) (?0) 


+ w. 
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ORIGINAL PAGE IS 
OF POOR QUALITY 

The parameters & N Q and represent the net change in dN°/d£ and 
dY 0 i 3 /d^ across the reaction zones at and respectively, and 

p(£) is the Favre probability density of £. In deriving these 

equations, the usual high Reynolds number assumptions have been 
made and the scalar dissipation is assumed not to be correlated 
with the mixture fraction. 

Finally, the mixture fraction pertubation, 

€* = € - I 


gives the transport equation for the variance of the mixture frac- 
tion 



9x k ' 


d ( £eff ag» 2 \ 

3X k Y ~2 dx k J 






( 91 ) 


The main features of measured probability density functions are the 
strong spike associated with the free stream and the continuous 
distribution generated by turbulence. By splitting the two, a 
clipped Gaussian intermittent formula is used to represent the 
Favre averaged probability density at a particular value, £ 

a 



y~ I 


for <0.25 t L 




1.25 




for £” 2 ^ 0.25 T 2 
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and 


original 9?m m 

OF POOR QUALI T1 


(? 2 ) = y [(? 2 ) f + (? t -f) 2 ] + (, -r)f 

( 93 ) 

Also from simple fits to the partial equilibrium calculations, the 
mean density as a function of the mole number and fuel mass func- 
tion pertubations is 



where v° ( £ ) is the Favre average specific volume for zero per- 
turbations. 

Lastly, the kinetic dependent source erms w” and w“ are given by, 

•* JT 

"w n = T JJJ q (£, n, y) (£, n, y) ( I /p) p (£ , n, y) d £ d n d y 

(95) 


and 



P JJJs* (£, n, y) p (£, n, y) d y d n d£ 


where 



- 9000 n 


2 moles 
kg sec 


( 96 ) 


from consideration of the partial equilibrium calculations and 
q(£, n, y) is a quenching function to allow for breakdown of the con- 
strained equilibrium represented as follows: 
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for £ < 0.01 


q (£> n, y) = B*(£) = 0 

=; 50 (£ - 0.0 1 ) 0.0 i < i* 0.03 

s I ? > 0.03 


Simplification of the joint P.D.F. yields 

= - p B*9000n 2 


The fuel source term w is evaluated based on the kinetic rate of 

29 y 

Duterque over the range of expected values of £, n, and y and again 
the joint P.D.F. is avoided, as 


where 



P 


\ 

p*«,n,7) P (#)df 


o 



8.53.10 


1 1 \( Y \ 

I 2 Jexp 



sec 


(97) 


where X02 oxygen mole fraction and temperature, T, are taken from 
the partial equilibrium calculations for Methane. 


The above three differential equations are combined with the 
other five Favre averaged equations, as mentioned above, written in 
cylindrical coordinates and transformed using the stream function 


3 ^= p'urdr 

and put in finite difference form using a Crank-Nicholson central 
difference .scheme. The nonlinear coefficient terms and source 
terms were evaluated as the mean of the upstream value and the 
first estimate of the downstream value. 

Then, for any "output" position in the flame, the effect of 
pertubations on species mole fractions and temperature is calculated 
from the local P.D.F., 
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I 

x. - f X. «,n,p) P 

o 
I 

T = J T (f.n.y) P <f) ^ 

O 

again, where joint probabilities have been avoided. 

3.4 Spray Evaporation/Combustion Models 

Since the influence of the liquid fuel spray on combustor per- 
formance is quite pronounced, an accurate spray model is essential 
for any combustor modeling effort. The modeling of liquid fuel 
sprays is discussed in this section. 

The spray model currently used by Garrett is based on a 
Lagrangian discrete-droplet approach allowing for droplet slip but 
no turbulent dispersion, Eulerian (continuous formulations) ver- 
sions allowing for dispersion, with or without droplet slip, have 
also been developed by Garrett. Both approaches offer advantages 
in certain circumstances. 

3.4.1 A Review of Spray Models 

A number of spray evaporation/combustion models have been dev- 
eloped and used with varying degrees of success. Several review 
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papers are available in the literature; ' an excellent recent 

32 

review paper is by Faeth. Faeth has divided the spray model work 
into two major categories: 

o Locally homogeneous flow (LKF) models 

,o Separated flow (SF) models 


(98) 

(99) 
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In LHP models, the gas and liquid phases are assumed to be in 
dynamic and thermodynamic equilibrium at all times, with no droplet 
slip. Consequently, the use of LHF models should be limited to 
finer droplet sprays. Although the LHF predictions tend to 
approach measurements as the droplet size reduces, the agreement is 
relatively poor, even for the sprays that have an SMD of around 30 
microns. Compared to the SF models, the LHF models are easy to 
use because they require minimum information concerning fuel injec- 
tor characteristics, fewer empirical constants, and shorter com- 
putation time. The LHF models give useful qualitative descriptions 
of the spray development, the rate of which is generally over- 
estimated. 

In SF models, finite interphase transport processes of mass, 
momentum, and heat are taken into account; and these models are 
therefore of more interest in gas turbine combustor modeling. The 
SF models can be broadly divided into the following two major cate- 
gories: 

o Discrete Droplet Models (DDM) - Lagrangian 
o Continuous Formulation Models (CFM) - Eulerian 

Both categories contain features that make their application to 

i 

practical combustors desirable. In DDM, the fuel droplets are as- 
sumed to exist sufficiently removed from each other that droplet- 
to-droplet interaction can be neglected. This assumption is quite 
reasonable for regions in the combustor other than very ''close to 
the fuel nozzle spray origin. Thus, the analytical modeling of a 
single droplet can be applied to the gas turbine spray that greatly 

simplifies the formulation. The nozzle spray is divided into a 

33 

number of size groups usually determined experimentally, with one 
droplet representing the behavior all the droplets in its group. 
A DDM is constructed for each size jroup. Given an initial velo- 
city and temperature as determined from the injector characteris- 
tics, the droplet trajectory is calculated through the flow domain 
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as governed by drag and other forces, until the droplet evaporates 
or exits the calculation grid. At each point along the flight path 
as evaporation occurs, modifications are made to the momentum, 
enthalpy, and species equations that govern the gas phase flow. 

The second major category of SF models is CFM. These models 
solve turbulent transport equations for the motion of the droplets 
and the turbulent diffusion of droplets is included. An underlying 
assumption is that all the droplets and the gas phase have the same 
velocity. The computational effort required for CFM is greatly in- 
creased because a complete equation (similar to the momentum, spe- 
cies, etc., equations) must be solved for each droplet group. Com- 
puter storage must also be allocated for the extra variables. A 
major disadvantage of this approach is that errors are generated in 
the vicinity of the fuel nozzle. Since the difference in the 
liquid and gas phase velocities is very significant in this region, 
a better resolution of grid spacing is needed than can be managed. 

3.4.2 Garrett's Spray Models 

The main requirement of a spray model is accurate predictive 
capability within a reasonable amount of computational effort, es- 
pecially for 3-D flows of practical interest. To achieve accuracy, 
various physical processes must be incorporated into the model in a 
realistic manner. Thus, relative velocity differences between the 
gaseous and liquid phases (droplet slip) , resulting in interphase 
momentum transfer, must be considered. Also, the evaporation of 
droplets during heat-up time (interphase heat and mass transfer) is 
important in order to predict ignition processes. Finally, tur- 
bulent diffusion of droplets is important but is generally ignored 

in most of the spray models. Stochastic models to consider tur- 
. 35 

bulent diffusion, as reported by Gosman are computationally ex- 
pensive when applied to real combustors. 
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On the computational aspect of spray modeling, it should be 
realized that a different degree of resolution is required in the 
near-injector and far-injector regions. 

Consideration of these factors has led to the development of 
the Garrett spray models, the features of which are described 
below. 

Evaporation during heat-up time is considered in the Garrett 
model and the computation of heat-up and evaporation rates includes 
realistic properties of jet aviation fuels. The spray model is ap- 
plicable to both dense and sparse sprays and is coupled into either 
the 2-step or the 4-step global hydrocarbon oxidation scheme; it is 
available in both 2-D and 3-D combustor performance pro- 
grams. 

Velocity differences (slip) between the droplets and the gas 
phase are modeled, and so is turbulent diffusion of the droplets.’ 
For computational purposes, the droplet sizro distribution is dis- 
cretized, _ differential equations in a Eulerian framework are 
solved for the velocity components and concentration of droplets in 
each size group (typically five size groups are considered) . To 
obtain good resolution in the near injector region, Lagrangian 
equations of motion for the droplets are solved in this region. 
Each class of droplets is tracked through the flow field in the 
vicinity of the injector and the interphase transports of mass, 
momentum, and energy are used to couple this solution to the 
Eulerian equations. The model thus combines the desirable features 
of the DDM and CFM approaches. 
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Lag rang i an Model 


The Lagrangian model was initially pursued since the formula- 
tion allows the tracking of small fuel droplets that are signifi- 
cantly smaller than typical grid dimensions and since the computa- 
tional requirements for 3-D flows are quite small. In addition, 
the assumption of sparse sprays and no droplet-droplet interaction 
is quite reasonable for most regions in a gas turbine combustor. 


The generalized governing equation for fuel mass fraction is 


— » III 

div <pVm fu - r fu eff grad m fu > = rh fi - R fu 


( 100 ) 


where R fu is the destruction rate of fuel and m e is the rate of 
vapor production from the fuel droplets. The vapor production rate 
or evaporation rate is determined from the burning rate constant 
k Q # which relates the change in the square of the droplet diameter 
(D) with time. 


k o - MUB) 


( 101 ) 


where: 


X-^ = Thermal conductivity of vapor 

Cp^ = Weighted average specific heat of vapor and air 
B - Mass transfer number 
- Fuel density 


From the burning rate constant, the fraction of fuel evapora- 
ted from a group of droplets (AF^,) can be determined from the time 
ihtegral; 


A Fj .= 


1.5 


P* D s 


r/ p f D r k i dt 


( 102 ) 
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where is the initial droplet diameter and the i subscript de- 
notes parameters of the ith droplet group. 


The total evaporation rate is then given by 


fu v 

m e = ~mr iF i 


(103) 


where: 


= The total fuel flow rate 
= Number of droplet groups 

= Volume of grid cell through which droplet is pass- 


The calculation of ft given above is based on the model of 
30 e 

Williams, which uses the temperature difference between the drop- 
let and ambient as the driving force. An alternate approach is 
that of Priem and Heidmann, who use the partial pressure of the 
fuel vapor as the driving force. The advantage of the Priem and 
Heidmann technique is its applicability to low temperature situa- 
tions sygfe altitude ignition, unlike the Williams model. 

^^pression for evaporation rate for the Priem and Heidmann 
model ■' ' W 


111 o 

rh = ttO l K Pa 
e vap 


(104) 


where: 


K = Function of vapor diffusivity and droplet Reynolds 
number 


* Vapor pressure of fuel at droplet surface 


vap 


a = 


In 


vap 


00 


P< ° “ P vap 


Poo - Ambient pressure 
Eulerian Model 


The salient features of the model are (a) velocity differences 
(slip) between the droplets and gas phase are allowed; and (b) tur- 
bulent diffusion of the droplets is included. For computational 
purposes, the droplets are considered to be present in a certain 
number (typically 5 to 10) of discrete size ranges. Differential 
equations in an Eulerian framework are solved for the velocity com- 
ponents and concentrations of the droplets in each size group. In- 
terphase transport of mass and energy (due to droplet evaporation 
and combustion) and momentum (due to drag between gas and liquid 
phases) are taken into account. Turbulent diffusion of drops is 
treated as if the droplets were present as a gaseous constituent. 
The approach used here is to assume the diffusion to be governed by 
a Fickian type law with an appropriate turbulent Schmidt number, 
assumed to be uniform over the flow field; but this is easily ex- 
tended by specifying the turbulent Schmidt number as a function of 
local flow characteristics. The model is applicable to both dense 
and sparse sprays insofar as the volume occupied by the droplets is 
included in the governing equations. 

The partial differential equations governing the droplet mo- 
tion and concentration of each size group are all written in one 
general form as 

-ft (R p<f>) + div {r (juft- grad 4>) } = S 0 (ig5) 
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where : 


R = Volume fraction of the droplet size considered 
P = Density of the droplets 

<f> = Velocity component or concentration of droplet size 
considered 

U = Velocity vector of droplet size considered 
M = Turbulent viscosity 
= Effective Schmidt number 
= Source of <£. 

The term S^ contains the pressure gradient, surface friction, 
and interphase drag if <f> is a velocity component; interphase mass 
transfer (evaporation rate) if is droplet concentration; and in- 
terphase heat transfer (heat-up) if <f> is the enthalphy. 

A drawback of the Eulerian model is that it cannot give ade- 
quate resolution in the near injector region* An excessively fine 
finite-difference grid would be required to obtain adequate resolu- 
tion. The Lagrangian method is capable of providing this resolu- 
tion in the near injector region. Garrett has therefore coupled 
two methods in order to utilize the advantages of each. The method 
is described next. 

Lagrangian/Eulerian Model 

The concept of this model is analogous to the near-wall treat- 
ment described in Paragraph 3.2 and it can be called the near-noz- 
zle spray treatment. The interaction between the Lagrangian and 
Eulerian solution is through the boundary conditions and the source 
terms . 

To obtain good resolution in the near-injector region, a spe- 
cial treatment is used in this region. The gas properties eval- 
uated by the Eulerian solution are used to solve a set of 
Lagrangian equations of motion for the droplets again allowing for 


if ~'x 
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interphase mass, momentum, and energy transfer. Each class of 
droplets is tracked through the flow field in the vicinity of the 
injector and the interphase mass, momentum, and energy transport 
that appear as source terms in the gas phase equations are thus 
evaluated. These terms are obtained by summing overall droplet 
size groups for each elementary control volume. The Eulerian equa- 
tions are then solved again including the interphase transport 
items. The procedure is repeated until it converges to the desired 
degree of accuracy. 

Physical Processes 

The physical processes involved in. spray modeling are inter- 
phase momentum transfer (drag forces) , interphase heat transfer 
(droplet heat-up) and interphase mass transfer (droplet evapora- 
tion) , and turbulent dispersion of droplets. A brief description 
is provided in the following paragraphs. 

Drag Forces - To calculate the drag forces on the droplet, the drag 
coefficient, C n , must be determined. Several expressions are 
available and Garrett has adopted the suggestion of Briffa who 
measured water droplet velocity decay using a shadowgraph tech- 
nique. Other forces, such as buoyancy or gravity, acting on the 
droplet, are quite small in comparison to drag and are usually 
neglected. 

Droplet Heat-Up and Evaporation - Calculation procedures for the 

rate of phase change of droplets fall into two basic categories: 

two^-stage and transient heat-up models. In two-stage models, as 

30 . 

discussed by Williams, the droplet is assumed to heat up to the 
boiling temperature with no evaporation occurring. Once obtained, 
the evaporation rate is governed by expression for the burning rate 
constant, defined as the time rate of change of the square of the 
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droplet diameter. Expressions can be used that account for the 
existence of an envelope or wake-type flame. The driving force is 
the temperature difference between the droplet and the surrounding 
gas phase. 

36 

Transit! t heat-up models such as that of Priem and Heidmann, 
use the difference in fuel vapor concentrations between the droplet 
surface and its surroundings as the driving force. The temperature 
of the drop is determined from the consideration of heat transfer 
to the drop and the fuel latent heat of vaporization. 

Though it offers the advantage of ease of computation, the 
two-stage models best represent droplets in the high-temperature 
zones of the combustor, where droplet heat-up time is quite short 
and local fuel concentration is low. 

Droplets exist for a significant period of time in the rela- 
tively cool, fuel-rich zone near the nozzle. The transient models 
better represent such droplets. The transient models are more com- 
plex from a computation standpoint, but they reflect the varying 
boiling temperature through the droplet life history and are sup- 
erior in predicting the evaporation in low-temperature environments 
(during an altitude start, for example) . j 

To evaluate the evaporation rates, fuel properties such as 
specific heat and distillation curves are required. For typical 
aviation fuels, these properties are usually available in the lit- 
erature or can be estimated from basic characteristics such as spe- 
38 

cific gravity. 

Turbulent Diffusion of Droplets - In most of the spray models, tur- 
bulent dispersion of droplets is ignored or introduced in an over- 
simplified manner. Some recent studies have adopted a stochastic 

35 

approach to model this feature. Recently, Gosman, et al., have 


68 


presented a stochastic discrete droplet method. In this method, a 
statistically significant number of random droplet samples is 
tracked in a Lagrangian framework and the ensemble-averaged 
behavior is assumed to represent the turbulent dispersion of drop- 
lets. This procedure is likely to be computationally expensive for 
real combustors where a large number of samples is required to 
obtain statistical averages. 

The Garrett Eulerian model includes turbulent diffusion. The 
model of diffusion of droplets is assumed to be the same as that of 
the gaseous phase; the extent of diffusion is controlled through 
the specification of the turbulent Schmidt number for droplet dif- 
fusion. 
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3.5 Soot Formation and Oxidation 


In this paragraph, soot formation and oxidation in combustion 
chambers are discussed. A general background on soot emissions is 
provided first. Quasi-global expressions for soot formation and 
oxidation a/..- described. A description of the influence of turbu- 
lence on soot formation and oxidation is included. The current 
approach adopted by Garrett is described next. This approach con- 
siders the influence of turbulent fluctuations on soot formation 
and oxidation rates. 

3.5.1 Background 

The particulate emission of primary concern in the combustion 
of hydrocarbon fuels is soot, which is evident in the form of 
exhaust smoke. The emission of smoke from gas turbine engines is 
responsible for the following problems of concern in this program: 

o Higher liner temperatures due to increased radiative heat 
transfer 

o Impingement of carbon on metal surfaces, resulting in 
erosion and reduced equipment lifetimes 

o Distortion of fuel spray distribution due to carbon 
deposits, leading to hot spots. 

Recently, attention is being directed toward the combustion of 
alternate fuels derived from coal liquids and shale oil. Since the 
use of these fuels results in significant increases in smoke pro- 
duction, a better understanding of the physical and chemical pro- 
cesses governing soot production is needed. 
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The processes governing the formation and subsequent oxida- 
tion of soot are of a particularly complex nature and, as such, 
quantitative models of soot production have yet to be developed. 
Soot is not an equilibrium product of combustion and, therefore, 
its formation is influenced as much by the physical processes of 
atomization, evaporation, and fuel/air mixing as by reaction 
kinetics. Soot is generally produced anywhere within the combustor 
where fuel/air mixing is inadequate, resulting in oxygen-deficient, 
high-temperature zones. 

For the pressures and temperatures normally prevalent in gas 
turbine combustors, equilibrium calculations indicate that solid 
carbon appears when there is insufficient oxygen to oxidize the 
hydrocarbon to CO and H 2 according to the relation: 

C H + * 0, - xCO + | (106) 

That is, the carbon-oxygen mass ratio for incipient soot formation 
is 12:16, or alternatively, the atomic C-0 ratio is unity. How- 
ever, since soot formation is essentialy a nonequilibrium phenome- 
non, experimentally, soot is observed at C-0 ratios (a) much less 

than unity at low temperatures (<2000 °K) ; and (b) greater than 

39 

unity at higher temperatures. 

Smoke levels are primarily dependent on 

o Air/fuel mixing 
o Temperature 
o Equivalence ratio 

o Residence time of air/fuel mixture 

o ' Pressure 
o Fuel composition. 
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These factors influence both the formation and subsequent oxidation of 
soot and are dependent on engine operating conditions, details of the 
combustor internal flow field, fuel droplet characteristics, etc. 


3.5.2 Quasi-Global Models of Soot Formation and Oxidation 

Since the elementary steps in the formation and oxidation of soot 

are not totally understood, Garrett uses quasi-global models that 

characterize soot production occurring via a few overall steps. Such 

40 

models have been successful in predicting soot production. 

The quasi-global models do not predict the size of soot par- 
ticles. W?th the current state-of-the-art, it is not possible to 
predict the size of formation of the soot particles in any practi- 
cal flow situation. Therefore, it is assumed that particles are 
produced at a known size in any analysis. It may also be assumed 
that particles are produced in accordance with a specified size 
distribution (e.g., Gaussian). 
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Tesner, et al., proposed a soot production model that 
grouped the complex processes of pyrolysis, nuclei formation, and 
soot formation into three rate-limited subglobal steps: 

Pyrolysis? 

n 0 = a o C fy ex P {parT./rn“.s) (107) 


Nuclei Formation: 
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Soot Formation: 


(108) 


(109) 


R sf = m D (a - bN)n (kg/m 3 .s) 


where a Q , E, g, g Q , a, and b are constants for given fuel; n Q is the 
rate of spontaneous formation of nuclei? n is the nucleus concen- 
tration; N is the concentration of soot particles? and m is the 

J» 

mass of a soot particle. 
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Khan and Greeves proposed a single-step global expression as 
a function of the partial pressure of unburned hydrocarbons (P HC )» 
the unburned equivalence ratio (0 u ) , and the temperature (T) : 

dC >3 o 

— = 0.468P|_j£<£!j exp (-40,000/RT) gm/crn s. (110) 


In both the above models, soot oxidation rates are not considered, 

Edelman, et al,^ consider both soot formation (R f ) and soot 
oxidation (R QX ) and express the net soot formation rate as 




-A f R 


ox 


( 111 ) 


where A fc equals total surface area available for oxidation. The 
formation step is expressed by a modified Arrhenius type of rela- 
tion: 


R f = AT a C|_^CQ^exp (-E/RT) gm/cm 3 s. (112) 

where Cc>2r C HC equal the concentration of unburned oxygen and 
hydrocarbon (gm/cm^) , and where A , a , a, b, E are model constants. 


C " 
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For the oxidation st*o, Edelman, et al., adopt the semi- 

43 

empirical formula of Nagle and Strickland-Constable for 
pyrolytic graphite oxidation; this formula is nonlinear and non- 
Arrhenius in P 02 an< ^ T: 


A t R OX - 12 


where 


/ k a p o \ 

[Wo f ^V-Wj 

A t gm/s 

(113) 

□ +k t / (k b p o 2 ):i ~ I 


(114) 

K a = 20 exp(-30,000/RT) 


(115) 

K b = 4.46 xI 0” 3 exp (-15,200/RT) 

(116) 

Kj= 1 .5 ! x 1 0 5 exp(~ 97,000/RT) 

(117) 

K z = 2 1 .3 exp(4 1 00/RT) 


(118) 


£ A i 

Shock-tube measurements of soot oxidation rates qualitatively 

confirm the features of the above formula. With these expressions 

for soot formation and oxidation and assuming a single soot part- 

0 40 

icle size of 250A* Edelman, et al., ' obtained close agreement of 
the predicted soot -concentration (mg/1) with the experimental data 
in a jet-stirred reactor. Thus, these expressions assume perfect 
mixing. In a gas-turbine combustor, however, regions of unmixed 
species will exist, and turbulence will also influence the soot 
production rates. As such, modifications to these f repressions are 
required before they can be used for a general 3-D turbulent flow. 
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3 C 5.3 Influence of Turbulence on Soot Formation and Oxidation 


45 46 

Magnussen, et al,, ~ ' have proposed a model that accounts for 
the influence of turbulent fluctuations on soot production rates. 
In turbulent flows, chemical reaction occurs when reactants at a 
sufficiently high temperature are mixed at the molecular level. 
The molecular mixing process is analogous to the dissipation (e) of 
turbulent kinetic energy k and is associated with the smallest 
scales of turbulence. Dissipation is concentrated in highly 
strained regions of the fluid occupied by fine structures with 
characteristic dimensions of the same magnitude as the Kolmogorov 
microscale. The reactants are molecularly mixed in these fine 
structures, where reaction occurs. Magnussen, et al., proposed the 
following expressions for the mass fraction contained in the fine 
structures: 


-3/4 

7* = 9.1 • (R f ) 


(113) 


where R fc is the turbulence Reynolds number, and the rate of 
transfer of mass per unit mass between the fine structures and the 
surrounding fluid is 

-1/4 

m=23.6-(R t ) £ (120) 

The rate of reaction is proportional to ft x where X is the 
fraction of small-structure eddies that are sufficiently heated to 
react. It is assumed that X is proportional to the ratio of local 
reacted fuel concentration and total fuel concentration. Thus, the 
rate of reaction is 


Rfu = 23.6 (R t r ,/4 « X C min (kg/m 3 .s) 


( 121 ) 
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where 



V<‘ + o 

V«l '+ T ) VC fu 


Cp f = Product concentration 
Ca, = Fuel concentration 


(122) 


C min ** the smailer of c f u (C02/i) and i is the stoichiometric 
oxygen requirement. The temperature T* of the reacting fine struc- 
tures is T above the local time-mean temperature T: 


where: 


AHp C in 

<if* u «r a t t f» mm 
T* = T + AT = T + ■ p ~£-~ — 

H r « the heat of reaction 
C - the specific heat. 

r 


(123) 


, and the surrounding temperature T° is 


T° = T- AT 


y*x \ 
I - y*xj' 


(124) 


Using Equations (107) and (109), the mean rates of nuclei and soot 
formation are then expressed as 


R n ,f = n o,T* * y * • X • P/P* + n 0jT o (I - y*-X) • P/P° + <f -g) n 

« 

-g Q n*N*y*x Pip*- g 0 n° N°(l -7*x) P/P° < 125 > 

R s , f = m p (a-bN^n^W// 5 * 

••+ m p . (a - b. N°)*n° (I Vy**> P/P° (126) 
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Finally, the mean rates of nuclei and soot oxidation are expressed 
as: 

R n,c = R f u n/C fu ( P art / m * s > (127) 

R s,c = R fu C s /C fu (k 9 /m3 * s > (128) 

Magnussen, et al., used this model to compute the soot concentra- 
tions in a turbulent C 2 H 2 diffusion flame. By adjusting the part- 
icle diameter [entered as m , the particle mass in Equation (109)], 

XT 

and the constant a Q in Equation (107) , good agreement with experi- 
mental measurements was obtained. 

3.5.4 The Garrett Soot-Emission Model 


The model adopted by Garrett for computing soot emissions 
under NASA Contract NAS3-22542 is described in the following para- 
graphs. 

The computation of soot emissions involves the solution of two 
additional transport equations for the concentrations of nuclei and 
soot. To complete the equation specifications, the source terms 
and the Schmidt numbers for these two variables are as follows: 

The source term for nuclei concentration is expressed as 

n,f n,c 

where R p is given by the smaller of the two values from Equations 
ri f iz 

(108) and (125), R_ _ is given by Equation (127). Thus, these 

n,c 

expressions amount to the use of the turbulent reaction rates, 
subject to the limitation that they cannot be greater than the 
rates under well-stirred reactor conditions. 
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The source term in the soot concentration Equation (111) is 
similarly expressed as: 

R s,f " R s,c 

where R_ - is given by the smaller of the two values from Equations 

5 f £ 

(112) and (126) ; R_ _ is given by the smaller of the two values from 

S f v 

Equations (113) and (128) . 

The turbulent Schmidt numbers cr and or for soot and nuclei 

s n 

concentrations are assumed the same as for gaseous fuel (i.e., 
0.9) . 

This model has been incorporated into the Garrett 3-D com- 
bustor performance program. Preliminary computations indicate its 
ability to make qualitative predictions. 
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3 . 6 Radiation Modeling 


An adequate treatment of radiative heat transfer from combus- 
tion products is essential for the prediction of gas-turbine liner 
temperatures and heat-transfer rates. For this purpose, Garrett is 

at present using the six-flux radiation model based on the 

47 

Schuster-Hamaker approximation. The influence of soot, CO 2 , 

and H 2 0 on the radiation properties (absorptivities and emissi- 
vities) is included in these equations. 

3.6.1 The Flux Methods 


In the flux methods, the angular distribution of radiation 
intensities is replaced by a number of discrete intensity vectors 
in different directions, thus reducing the complexity of the 
integro-differential equation of radiation heat transfer. The 
energy transfer in each direction is represented by a closed first- 
order ordinary differential equation obtained by integrating the 
radiation transfer equation over a solid angle. This method was 
originated for the 1-D case as the two-flux method, wherein only 
two directions are considered. Considerable errors exist in the 
two-flux solution in the case of essentially 1-D heat transfer bet- 
ween parallel plates; a situation for which the method is supposed- 
ly best^suited. This suggests that the two-flux method is not 
sufficiently accurate to permit its application to the prediction 
of radiant transfer in practical systems. 

48 

Spalding extended the 1-D formulation to tvro and three 
dimensions by formulating the four- and six-flux models. Exten- 
sions of the two-flux model to multi-flux and nongrey emitting 
absorbing media are also discussed by Siddall. The four-flux 
model applied to an axisymmefcric combustor underestimated wall 
radiation fluxes, although temperature predictions were reason- 
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The reasons for the inaccuracies in these flux methods are (a) 
the radiant flux is divided into too few directions (2, 4, or 6 
being small for many applications or (b) the fluxes in the differ- 
ent directions are unrealistically independent of each other. 
Another limitation of the flux models is that their extension to 
general curvilinear coordinates for handling complex geometries is 
rather cumbersome. 

3.6.1. 1 The Six-Flux Model Used at Garrett 

A six- flux radiation model based on the Schuster-Hamaker 
47 

approximation is used currently at Garrett. It should be noted 

49 

that, as pointed out by Siddall , other flux model approximations 
such as Milne-Eddington and Schuster-Schwarzschild can be repre- 
sented by the same form of flux equations with constants being 
different. 

The differential equations describing the variations of the 
fluxes along six directions can be reduced to the following three 
second-order ordinary differential equations: 


s; ( 5TT o(R x -E).|(2R x -R r -R Z ) (129a) 

r$ r( 3T!Tl-f- ) = a < Rl '- E ) + ! (2R r -R x -R Z ) (129b) 

r 

F3F <^r *7W > * 0 (rZ - e > ♦ t (2rZ - rX - Rr> (129o) 
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Where the composite-fluxes R , R and R are defined as 


R M «x + + 'x-> 

R r 4 ( l r+t l r .) 




where I x+ , I r+ and 1 Q+ are the fluxes along the positive directions 
of axial, radial, and circumferential directions, respectively; 
I x _, I r _, and I ^ are the corresponding fluxes along the negative 
directions. 


a = absorption coefficient, defined as radiation absorbed 
per unit length 

s = scattering coefficient, defined as radiation scattered 
per unit length 

E = black body emissive power = a 

cr = the Stefan-Boltzman constant. 

3.6.2 Discrete Transfer Method 
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Lockwood and Shah have presented a method called the 
discrete transfer method. This method is based on the solution of 
representatively directed beams of radiation within the combustor, 
as in the Monte Carlo method. However, in this method the direc- 
tions of the rays are specified in advance and they are solved for 
only between two boundary walls contrary to the Monte Carlo method 
where the ray directions are specified at random and the rays are 
tracked to extinction. Lockwood and Shah have shown that this 



method closely reproduces the analytical solution for radiation 
between parallel plates (1-D case) , radiation in a square enclosure 
(2-D), and in a cubic enclosure (3-D). The conventional two- , 
four- , and six-flux models for these cases show larger errors. 
This new method is economical, geometrically adaptable, provides 
ease of application, and has the possibility of obtaining any 
degree of precision (through the specification of number of rays). 
The method is designed to be coupled to fluid flow solutions. GTEC 
has used this method in its 2-D combustor program. 

3.6.3 Radiation Properties 

The contributors to radiation fluxes in gas turbine combustors 
are: soot, CG 2 , HjO (vapor), inorganic particles, etc. Only the 
influence of soot, C0 2 , and H 2° (vapor) is discussed here. 
Although CO and unburned C H r contribute to emission and attenua- 
tion of radiation within flames, these contributions are localized 
and of secondary importance for computing radiative fluxes. The 
contributions of N0 x and S0 2 can be neglected because of their low 
concentrations . 

The radiation properties of the principal radiating species 
including soot, C0 2 , and H 2 Q are significantly nongrey. Conse- 
quently, the calculation of the radiation properties > is a time- 
consuming task. However, detailed spectral calculations are 
unnecessary since approximate calculations (by means of curye fits) 
are more convenient and provide good accuracy. ^ Garrett has 
employed the approximate curve-fit procedure for the calculation of 
radiation properties under NASA Contract NAS3-22542. 14 

The absorptivity (a) of the gas-soot mixture includes the soot 
absorptivity, the absorptivity due to the absorption bands of C0 2 , 
and H 2 0 and corrections for the overlapping of bands. 


82 


53 


Using the spectral data,*''*' the gas absorptivity is calculated 
by taking a summation over the absorption bar*ds of C0 2 and H 2 0. In 
the approximate calculation method adopted by Garrett# a simpler 


approach is used. 


The gas absorptivity# oig is written as 
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a 


9 e g 


(T/T ) 


where £ = P W /(P W + P c ) 

e = gas emissivity at a temperature T and path length LT /T 
y s i 

\ 

f 

\ 

P c' p W = P art i a -*- pressure of C0 2 and H 2 0 I 

e g € c + e w ^ e cw 

where 

e c' e w = s missi vities of C0 2 and H 2 0 

Ae cw = overlap correction factor. 

e can be computed using a temperature adjusted version of 
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Leckner's approximate overlap correction # and approximating 

cw 

e and g by curve fits of P , p , T and path length to spectral 
c w c w 

calculations. In the range of interest in gas-turbine combustors, 
such calculations agree to within 5 percent of the spectral calcu- 
lations and the experimental results. The absorptivity («) of the 
gas-soot mixture is given by 


With a _ obtained above, it remains to determine a , the soot 
9 s 56 

absorptivity. This is obtained by the method of Felske and Tien. 

This method assists that the complex refractive index of soot is 

independent of wavelength and that the soot particle diameter is 

small compared to the wavelength of radiation, so that scattering 

is negligible. The spectrally integrated absorptivity can then be 

written in a closed-form expression to determine <* . 

s 

By using the radiative property calculations of the type 
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described above, Sarofim indicated that radiation calculations 
can be made with fair confidence, and that the major source of 
uncertainty in such calculations is soot concentration, rather than 
gas-radiatoh properties. 


) 
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SECTION IV 

4.0 DESCRIPTION OP THE COMPUTATIONAL SCHEME 

The governing differential equations described in Section 3.Q 
are nonlinear and coupled partial differential equations. In most 
practical situations, it is not possible to obtain analytical solu- 
tions to these equations and numerical methods have to be used. A 
description of the numerical scheme used in the CPM is provided in 
Paragraph 4.1. The treatment of the boundary conditions is given 
in Paragraph 4.2, and the criteria for convergence and the method 
for assessing grid independence are outlined in Paragraph 4.2. 

4.1 Description of the Numerical Method 

The numerical method used in the CPM's are based upon the 
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finite difference technique of Patankar, which used the Semi- 
Implicit Method for Pressure Linked Equation (SIMPLE) algorithm. 
The features of this computational procedure include the following: 

& Solution of a sufficiently general single form of differ- 
ential equations 

o Provision for use with different physical models 

o Use of pressure and velocities as the main hydrodynamic 
variables 

o Use of the pressure-correction technique 
o Use of nonuniformly spaced grids 
o Use of staggered storage locations 
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o Derivation of finite-difference equations by integrating 
the differential equations over finite control volumes 

o Ijine-by-Xine solution of the difference equations 

The finite-difference equations are derived for a box-shaped 
flow domain. Over the region of interest, a number of grid planes 
parallel to the two coordinates are placed. For each grid node, 
the finite-difference equations are set up for each of the flow 
variables to be solved. Since the governing equations for axial- 
and radial-velocities (Equation 1) contain pressure gradient terms, 
these two variables are solved along planes staggered with respect 
to the main grid planes described above. 

A typical grid node spacing for a general flow problem is 
shown in Figure 4.1-1. Finite-difference equations for a node are 
obtained by integrating the differential equations over a control 
volume enclosing a grid node. For evaluating the convection and 
diffusion fluxes through a control volume face, a linear variation 
(in the direction normal to the face) of the flow properties is 
assumed. For other purposes, a stepwise variation with discontinu- 
ities at control-volume boundaries is assumed. Net rate of flow of 
$ into the control volume around a node P (Figure 4.1-1) by convec- 
tion and diffusion in the x-direction is 

f T X- + ^ ■ f x-* L x- 3 < ^x- +[ T x+ " f x+ L X+^X+ 

‘ [T X-' f X- L X- + T X+ + a “W L X+ 3 ^P (130) 

where T x = T eff ^ x / 8X 
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m x ‘'/SX 


L X = 

A x = 0.5(r + + r_)AY 


V Defining/// S^dV = S u + S p tf p , 
equation for the variable becomes 


the one-dimensional transport 


[ T x _ + (1 - f x _) L 
= [T x _ + (1 - f x J L y _ 


X-t T x+“ f X+ L X+" S p ] 
^x- +[T x + ’ f x + L x + ^x + 


+ s 


u 


( 131) 

The linear-profile assumption becomes unacceptable when f x+ L x+ is 
large compared with T x+ because with weighting factor (T x+ - f x+ 
L x+ ) then becomes negative, implying an unrealistic physical pro- 
cess through which raising the value of </> x+ could lower the value 
of Therefore, it is assumed that if the convective flow rates 

(L) are large compared to the diffusion coefficients (T) , the dif- 
fusion across the control-volume face is zero and the value of 
convected is equal to the value at the node on the upwind side of 
the face. With this assumption, the coefficient T x+ - f x+ l x+ is 
replaced by T* + - F x+ L x+ where 

T X4 = tT X+’ ■ (1 - W L x+’ J x+ L x+' 

Here [a^, a 2 , a 3 ] stands for the largest of the three quantities 
a^, and a 3 * 


87 


Using a similar procedure for the fluxes in the radial direc- 
tion, the final finite-difference equation is reduced to 

A p*p = A X+X+ + A X- X- + a y-^y+ + A Y-^Y- + S u 

(132) 

The solution of the above equation is obtained by line-by-line 
relaxation using an efficient tri-diagonal matrix algorithm. By 
this method, a traverse along one direction, for example, the X- 
direction, is made with old values for the y-direction nodes. 
Using this solution as the best estimate, the y-direction is then 
traversed. The solution method adopted is based on the SIMPLE 
algorithm of Patankar and Spalding as described in Reference 15. 
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4.2 Boundary Conditions 


The boundary conditions are enforced by appropriately modify- 
ing the finite-difference coefficients at the first interior point 
adjacent to the boundary. For the inlet boundaries, the velocity 
components, density, and turbulence profiles are either experimen- 
tally known or estimated. At the inlet boundary, if pressure is 
specified, the pressure correction is set to zero. When the normal 
velocities at the boundary point are specified, the coefficients in 
the pressure correction are modified in such a way that the mass 
fluxes through the control volume satisfy the overall continuity 
equation. 

For boundaries of the second kind, where gradients and not the 
values of the variables are specified, the program uses one of the 
following two approaches* In the first approach, the boundary 
value is guessed and continually updated to satisfy the given gra- 
dient condition. The second approach breaks the link through the 
boundary to all adjoining external control volumes by first arrang- 
ing for the finite-difference coefficient connecting the boundary 
node to an internal node to be zero, and then inserting the correct 
flux at the boundary as a source of diffusion and/or convection for 
that internal node. 

At the symmetry plane, the convection and diffusion fluxes in 
the radial direction are zero. Therefore, the finite-difference 
coefficients containing these fluxes are set to zero at the axis of 
symmetry. For the exit plane, information about some of the vari- 
ables is not available. However, since it is the process occurring 
in the calculation domain that decides values of the variables that 
the outgoing fluid will carry, there is no need for information at 
such boundaries. These boundaries are simply treated by neglecting 
the diffusion at the exit boundary. 
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Boundary conditions at the near-wall nodes are treated in the 
manner outlined in Section 3.0 (Equations 7 through 9) . 


The input parameters depend upon the nature of flow problem 
computed. In many of the test cases, initial profiles of turbu- 
lence kinetic energy (k) and length scales (L) are not available. 
For these cases, uniform profiles of k and L are prescribed at the 
inlet and the default values used are 

k = 0.003 U?„ 


L 


0.02 


max 


where, U av is the average inlet velocity, and R max is the maximum 
cross-stream dimension of the flow geometry. If information about 
turbulence intensity levels is available at the inlet, appropriate 
uniform k values are used at the inlet. 
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4.3 Convergence Criteria 


The solution is accepted as the converged solution when the 
total mass source error is less than about 0.1 percent of the total 
mass flow rate. For all the test cases , the computations were car- 
ried out further to ensure that the profiles of the dependent vari- 
ables did not appreciably change. For all the test cases con- 
sidered in this program, when the solution converged to the accept- 
able limit of 0.001 on the total mass source error, the maximum 
mass source error in the computational domain was less than 0.0002. 
The number of iterations required to reach the acceptable conver- 
gence level varied from problem to problem. In most of the recir- 
culating flows, a minimum of 350 iterations were needed to reach 
the convergence criterion. 

During the computations, the values of each dependent variable 
are monitored to ensure that the maximum change in the value of 
each dependent variable is a small fraction of the reference value. 
When this condition is satisfied, and if the total mass source 
error is less than 0.1 percent, plots of all the variables of 
interest are obtained. Computations are then continued for another 
50 iterations and plots are obtained again. If these plots are 
identical to within graphical accuracy, the solutions are accepted 
as converged solutions. > 

The numerical solution obtained for any given flow problem 

j 

depends upon the grid density and grid distributions. The solu- 
tions are accepted as grid independent if the predicted results are 
essentially invariant when the grid density or the grid distribu- 
tions are changed. This type of test was performed for many of the 
test cases, but these test results will be presented only for a few 
of them. For the other cases, the predictions presented in this 
report are essentially grid independent. The details about the 
grid distributions for each test case will be provided along with 
the discussion of the results. 
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SECTION V 


5.0 DATA BASE FOR BENCHMARK QUALITY TEST CASES 

To assess and critique the current models and generate a pro- 
gram plan to improve their accuracy and usefulness as a combustor 
design tool, the assessment of the models was conducted, in the fol- 
lowing two interrelated steps: 

o Assess and critique the physical submodels involving the 
fundamental processes of combustion, individually, with 
data from ideal element tests under well-defined condi- 
tions. The physical submodels considered here are turbu- 
lence modeling, gaseous fuel combustion, spray evapora- 
tion and combustion, soot formation and oxidation, and 
radiation modeling. 

o Assess and critique the model predictions against the 
data from advanced gas turbine combustors. 

Accordingly, the data base is arranged in two sections: Para- 
graph 5.1 includes a description of the data base from ideal 
element tests and Paragraph 5.2 contains a description of the data 
from a number of gas turbine combustors. 

5.1 Data Base from Ideal Element Tests 

A literature survey of recently published work (generally 1970 
or later) was conducted to compile a data base necessary for a 
benchmark quality test case. Published literature related to the 
following submodels was reviewed: if \ 
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o 


Turbulence Modeling 



o Gaseous Fuel Combustion 
o Spray Evaporation and Combustion 
o Soot Formation and Oxidation 

In the following paragraphs, a data base from ideal element 
tests is provided. The ideal tests range from simple entrance 
flows in pipes and 2-D channels to more complex flows like the flow 
fields behind steps, blockages, and swirling recirculating flows. 
These tests are intended to encompass the range of complexities 
involved in combustor internal flows. Simple entrance flows are 
included in the validation efforts to ensure that the analytical 
models can be used to predict simple flows without any modification 
to the model. The data base selected has fairly detailed measure- 
ments, including turbulence parameter measurements, with estima- 
tions on errors. 

5.1.1 Turbulence Modeling 

In this paragraph, a data base for assessing turbulence models 
is provided. The assessment procedure for the k-e turbulence 
models will consist of comparing the predicted time-mean velocity 
components with the corresponding measurements. For the algebraic 
and full Reynolds stress models, the predictions of the turbulence 
intensities and cross correlations will also be compared with the 
measurements. Some cases involving scalar transport are also con- 
sidered, and these involve predictions of the concentration of a 
trace gas (inert) or temperature under heated but inert conditions. 

The references reported in the following tables provide infor- 
mation about the available measurements reported in the literature. 
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These references were selected based upon the extent and accuracy 
of the data and the nature of the geometry of test conditions. The 
references selected are presented in the form of increasing order 
of complexity of the flow field in the form: 

o Simple Flows (Boundary Layers, Jets, Mixing Layers, etc.) 
(2-D Parabolic) - Table 3 

o Streamline Curvature Effects (Curved Ducts, Curved 

Boundary Layers, etc.) (2-D Parabolic) - Table 4 

o Recirculating Flows (Nonswirling) (Both Un^nfined and 
Confined) (2-D Elliptic) - Table 5 

o Swirling Flows (With and Without Recirculation) (2-D 

Elliptic/2-D Parabolic) - Table 6 

o Scalar Transport - Table 7. 

From the references provided in these tables, benchmark test 
cases were selected, as described in Sections 6.0, 7.0 and 8.0. 
These cases were used to evaluate the turbulence and kinetic model 
predictions. 
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5.1.2 Gaseous Fuel Combustion 


In this paragraph, a data base for assessing gaseous fuel com- 
bustion models is provided. The assessment procedure consists of 
comparing the predictions of time-mean velocity components, temper- 
ature, concentrations of species {unburned fuel, CO, C0 2 , H 2 , H 2 0, 
0 2 , n 2 ) against the experimentally measured values of these quanti- 
ties. These quantities were selected because they are of interest 
in gas turbine combustors. Reliable measurements of these quanti- 
ties are available, and they are a good indication of the predic- 
tive capability of the gaseous combustion model consisting of the 
turbulence/chemistry interactions and the hydrocarbon reaction 
mechanisms. The assessment will be done for different flow types: 
turbulent/laminar , premixed/diffusion, one/ two/ three-dimensional 
flow, parabolic/elliptic, swirling/nonswirling. 

In accordance with the assessment procedure, the data base is 
categorized into four sections: 

o Laminar Premixed Flames - Table 8 

o Laminar Diffusion Flames - Table 9 

o Turbulent Premixed Flames - Table 10 

o Turbulent Diffusion Flames - Table 11. 

In each of these tables, the data is arranged in order of increas- 
ing complexity, starting from 1-D parabolic to 3-D swirling 
elliptic flows. 
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During the search for compiling the data base, several publi- 
cations were encountered wherein the boundary conditions or other 
information required for modeling were not clearly or completely 
stated. Such cases (e.g., References 111-123) have not been 
included here. Measurements of quantities not related to the 
assessment procedure given above have also been excluded. The data 
base is concerned with the measurements of quantities listed above 
for steady gaseous hydrocarbon flames. 


> 
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TABLE 11. TURBULENT DIFFUSION FLAMES 
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5.1.3 Spray Evaporation and Combustion 


In this paragraph, a data base for assessing spray evaporation 
and combustion models is provided. The assessment procedure will 
consist of comparing the predicted spray trajectory, droplet con- 
centrations, velocities and size distribution, temperature, and 
concentrations of species (unburnt fuel, CO, CC> 2 H 2 , H 2 0, 0 2 , N 2 ) 
against the experimentally measured values of these quantities. 
The available data on spray evaporation and combustion is listed in 
Table 12. The predictions of these quantities is an indication of 
the accuracy of the various features of the spray model: 

o The prediction of spray trajectory, droplet concentra- 
tions, velocities, and size distribution under non- 
burning and nonevaporating conditions reflects on the 
accuracy of the spray dynamics model, which includes the 
modeling of the drag forces between the spray and the gas 
phase. 

o The prediction of the droplet concentrations and size 
distributions along with the mixture fraction under non- 
burning (but evaporating) conditions serves to test the 
droplet heat-up and evaporation models. 

/ 

o Finally, the prediction of droplet concentration and size 
distribution along with gas temperature and composition 
serves to test the validity of the spray combustion 
model . 

Thus by assessing the predictions of the quantities listed 
above, all features of spray evaporation and combustion involving 
interphase momentum (spray dynamics, drag), feat (droplet heat-up) 
and mass (droplet evaporation and combustion) transfer are tested 
individually and jointly. 
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'fij.e data base is arranged in order of increasing flow complex- 
ity. Sources in which all boundary and initial conditions required 
for modeling were not completely or clearly stated have not been 
included in the data base (e.g. References 156-170) . It should be 
noted that complex two-phase slip models as used at Garrett require 
detailed information specifying the initial conditions at the fuel 
injector? initial drop size distribution, initial velocity distri- 
bution, etc. This information is generally not available and 
therefore has to be estimated from the available injector charac- 
teristics. 
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5.1.4 Soot Formation and Oxidation 


In this paragraph, a data base for assessing soot formation 
and oxidation models is provided. The assessment procedure will 
consist of comparing the predicted soot concentration, temperature, 
concentrations of species (unburnt fuel, CO, C0 2 , H 2 , H 2 0, °2' N 2^ 
against the experimentally measured values of these quantities. 
The comparison of predicted and measured soot concentrations is a 
direct indication of the accuracy of the soot model. Temperature 
and gas composition are affected by the presence of soot to an 
extent depending on its concentration. Therefore, assessing the 
accuracy of the predictions of temperature and gas composition 
serves to indirectly assess the ‘soot model. 

The data base for the soot models is rather inadequate since 
very few measurements under controlled conditions have been 
reported in the literature. The reason is the difficulty in accur- 
ately measuring soot concentration profiles in a combustor. Quite 
often, only the exhaust smoke concentration is measured and soot 
profiles have been measured in only simple flames. 

As in the preceding sections, several sources of data (e.g. 
References 185-194) were found that were not suitable for model 
assessment due to incomplete specification of the boundary and 
initial conditions. These have not been included here. Also, mea- 
surements related only to gas turbine type fuels have been con- 
sidered, since it is practically impossible to validate the model 
and obtain a set of model constants for all types of hydrocarbons. 
The data base for soot formation and oxidation is presented in 
Table 13. 
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TABLE 13. SOOT DATA (CONTD) . 
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5.2 Data Base from Garrett Gas Turbine Combustors 


A number of gas turbine combustors have been mapped at Garrett 
over the last ten ytears . A brief description of the Garrett data 
base is given in the following paragraphs. 

5.2.1 Can Combustor Mapping 

A nonreacting can combustor with swirlers at the dome was map- 
ped ' in 1973 for comparison with the k-e turbulence model. A 
schematic of the burner along with the flow split is shown in Fig- 
ure 5.2-1. A calibrated three-hole wedge probe and liquid micro- 
manometer were used to measure the radial distirubtion of the yaw 
angle, static and total pressures at different axial stations. 

As part of model validation under the USARTL Design Criteria 
Program, another can combustor nonreact,ing flow was mapped at dif- 
ferent throughflow rates. This combustor was filled with 21 mea- 
surement ports. 

A calibrated five-hole pyramid probe was traversed across the 
can combustor at three circumferential locations and seven axial 
stations. Four traverses were made in the primary zone, seven in 
the intermediate, and ten in the dilution zone. The probe mouncs 
and the test conditions are shown in Figure 5.2-2. 

Reacting flow mapping was accomplished on a similar can com- 
bustor, shown in Figure 5.2-3. Radial profiles of CO, CO.,, NO , 
and unburned hydrocarbons were measured at axial stations 6.0, 8.5, 
10.4, 12.9, 15.4, 18.8, 21.3, and 26.2 cm downstream from the fuel 
nozzle face. Five circumferential stations were mapped to deter- 
mine the profile variations in the circumferential direction. The 
mapping was conducted for both gaseous (natural gas) and liquid 
fuels (Jet A) over a wide range of operating conditions. The fuel 
nozzles used for each fuel ase shown in Figure 5.2-4. 
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REACTING FLOW CAN COMBUSTOR GEOMETRICAL DETAILS 

NO. OF GEOMETRIC AXIAL DISTANCE 


ORIFICE TYPE 

ORIFICES 

SIZE (cm) 

AREA, cm2 

(cm) 

DOME LOUVERS 

30 

0.36 

3.02 

_ 

PRIMARY 

6 

1.12 

5.89 

9.09 

DILUTION 

6 

1.42 

9.53 

17.21 

COOLING SLOT LIP 
#1 

30 

0.44 

4.6 

5.05 

#2 

30 

0.46 

5.43 

12.20 

#3 

30 

0.48 

5.43 

20.59 

#4 

30 

0.48 

5.43 

29.67 


Figure 5.2-3. Can Combustor for Reacting Flow Mapping. 
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Figure 5.2-4. Natural Gas Nozzle and Airblast Nozzle Used for 
the Can Combustor Mapping. 
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A water/steam-cooled stainless steel emissions probe with ten 
individual radial sampling points was used for the combustion emis- 
sions mapping. 

5.2.2 Annular Combustor Mapping 

A reverse-flow premix/prevaporizing (PM/PV) annular combus- 
tion system that is compatible with the Garrett TFE731-2 turbofan 
engine envelope was designed and tested to demonstrate combustor 
technology capable of meeting the 1979 EPA emission standards for 
TI class engines as part of the NASA Pollution Reduction Technology 
Program. To better understand the performance characteristics of 
this combustion system, internal radial profiles of gaseous emis- 
sions were measured in an atmospheric test rig. 

The piloted PM/PV combustion system incorporates two axially 
staged burning zones, as shown in Figure 5.2-5. The radial pro- 
files of CO, C0 3 , UHC, and N0 x were measured at four different 
axial-stations and six circumferential (0) planes within the main 
combustion zone. A water/steam-cooled probe was used to obtain 
radial profiles. The internal emissions mapping was conducted at 
one atmosphere in a combustor rig without the transition liner. 
The effect of different parameters including combustor inlet tem- 
perature (T 3 ) , overall fuel/air ratio, and fuel-flow splits between 
the pilot and PM/PV combustion zones on the emissions profiles were 
studied. The mapping was conducted with propane as PM/PV fuel to 
simulate complete evaporation; however. Jet A fuel was used for the 
pilot. 

Information concerning the internal flow field of a TFE731 
production combustor (Figure 5.2-6) was provided through measure- 
ments of CO 2 , CO, UHC, and N0 X taken inside the combustor primary, 
intermediate, and dilution zones at atmoshperic test conditions. 
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Figure 5.2-5. Axially Staged Burning zones of the Piloted 
PM/PV Combustion System. 
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The emissions probe used for the internal mapping was the same 
as the one used earlier to map the PM/PV Concept 3 combustor. The 
eight individual sampling ports of the probe were manifolded 
together to obtain only averages in the radial direction. The mea- 
surements were taken at different axial stations (as shown in Fig- 
ure 5.2-6) in the primary and secondary regions. 

Detailed internal gaseous emissions and temperature measure- 
ments inside the Uprate T76 combustor primary zone have been con- 
ducted at various axial locations. 

A single-point, water-cooled, emissions probe with an end cap 
(Figure 5.2-7) was designed for use in the primary zone. This 
probe is intended to separate relatively large liquid fuel droplets 
from the gas sample. The end-cap feature was also used in the con- 
struction of a ceramic radiation shield for the aspirated thermo- 
couple used to measure primary zone temperatures. 

Emissions samples were taken at five different axial positions 
from 1.016 to 6.35 cm from the dome. Temperature measurements were 
taken from the dome to the dilution zone. The measured sector 
extended over 36 degrees and was centered on a main fuel nozzle. 
Seven circumferential stations were selected to correspond with 
areas of carbon deposition in the Uprate T76 combustor. 

Most of the data was taken at an altitude idle-engine condi- 
tion and also at the sea-level design condition. Two combustors of 
the same part number were measured. Fuel/air ratios nearly twice 
the stoichiometric values were measured at the discharge of the 
primary zone for the design condition, indicating poor primary zone 
mixing . 
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Figure 5.2-7. Schematic of Emissions Probe and Measurement 
Locations for Fuel/Air Rate Profiles in the 
UT76 Combustor. 
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SECTION VI 


6.0 SIMPLE FLOWS 

The results for the benchmark test cases (shown in Table 3) 
are presented in the following four categories: 

o Model evaluation for simple flows 
o Model evaluation for complex nonswirling flows 
o Model evaluation for swirling flows 

o 3-D jet-mixing flow validation 

These categories are selected in increasing order of complexity 
and, for each category, the results will be presented for nonreact- 
ing and reeacting flows. In this section, discussion of results 
and model evaluation for simple flows are presented. To present the 
results, the predictions will be shown by lines and the data will 
be represented by symbols throughout this report. 

6.1 Flow Over a Flat Plate 


One of the benchmark test cases selected from the assembled 
data base is the flow over a flat plate, for which measurements 
were made by Watts and Brundrett . Their test plate was 2.44 m 
long with boundary layer trips placed near the leading edge to make 
the boundary layer fully turbulent. The mean velocity and the tur- 
bulence velocity fluctuations were measured with a hot-wire probe 
at x = 0.244, 0.462, 0.8466, 1.163, 1.4656 and 2.2743 m. The free 
stream velocity for this test case was 20.8 m/s. A schematic of 
this flow geometry is shown in Figure 6.1-1. 

Computations for this case were made using the Garrett 2-D 
parabolic code, and predictions were obtained with the following 
models: 
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o Standard k- e model 

o Standard k~e model with near-wall low Reynolds number 
correction 

o Algebraic stress model (ASM) 

o ASM with low Reynolds number correction 

For all these cases, the initial conditions were applied at x 
= 0.244 using the measured profiles. One hundred cross-stream grid 
points were used in these computations. The grid distributions 
were selected so that the nodes were closely spaced near the wall 
and are farther apart near the edge of the boundary layer. For the 
standard k-e model, the wall function treatment outlined in Section 
4.2 was used to specify the wall boundary conditions. 

The predicted mean velocity profiles using the standard k-e 
model are shown in Figure 6.1-2. This figure shows that the agree- 
ment between data and predictions was poor. Problems in this com- 
putation were associated with the wall function approach for pre- 
scribing the boundary conditions at the near-wall nodes. 

i 

One way to circumvent the application of the wall functions is 

to apply low Reynolds number corrections to the k and e equations 

that will enable k and e to be zero at the wall in a consistent 

manner. From the survey of literature for low Reynolds number cor- 

17 

rections, the model of Chien was selected for these computations. 
In Chien’ s model, the source terms and exchange coefficients in the 
k and e equations have been modified. The governing equations for 
k and € , still retain the form shown in Equation (1) . The differ- 
ence arises in equations (2) and (4) . The modified terms in 
Chien 's model are 
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u* is determined by using the linear part of the law of the 


u = u* y +J 0^ y + < 11.5 (140) 

With the modified terms in the k and € equations, the wall boundary 
conditions at the wall (y=0) are: 


k = 0 


e = 0 
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It should be recognized that for improving near-wall solution 
accuracy with Chien's mod if ication, one must employ a number of 
grid points inside the viscous sublayer. This is not always pos- 
sible for elliptic flows. Consequently, in this report, Chien's 
corrections are applied only for parabolic flows. In all the com- 
putations with Chien’s correction, approximately 10 nodes were dis- 
tributed in the viscous sublayer. 

A comparison between the data and the predicted results using 
Chien's correction are shown in Figure 6.1-3. A significant im- 
provement in the agreement and predictions is obtained with the low 
Reynolds number correction over the wall function approach in the 
standard k- e model. 

Figure 6.1-4 illustrates the predicted results for mean veloc- 
ity using the algebraic Reynolds stress model. In the ASM the 
value of the coefficient C D is computed from equation (23) while 
the standard k-e model assumes a constant value of C D . This param- 
eter, C D , is used for calculating the turbulence diffusion rate. 
Figure 6.1-4 shows that the ASM is in excellent agreement with the 
data of Watts and Brundrett. Application of the low Reynolds 
number correction on the ASM does not appreciably improve the mean 

velocity predictions. These results are presented in Figure 6.1-5. 

v 

However, the application of Chien's low Reynolds number correction 
substantially improves the prediction of turbulence kinetic energy 
components. A comparison of the predicted turbulence kinetic 
energy profiles at x « 1.8735 m using the four models are presented 
in Figure 6.1-6. The predicted turbulence kinetic energy values, 
using the standard k - e model, are slightly higher than the measured 
values. The peak k value near the wall is significantly higher 
than the measurements. When Chien's low Reyno.lds number correction 
is applied to the k- € model, the near-wall kinetic energy values 
are in better agreement with the data. The ASM predictions for k 
are also higher than the data, but is slightly better in comparison 
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with the k-€ model predictions. The application of the low 
Reynolds number correction significantly improves the predicted k 

values . 

A comparison between data and predicted time mean fluctuating 
velocity components at x = 1.8735 is shown in Figure 6.1-7. The ASM 
predicts much higher peak values for u', v' r and w' components com- 
pared to the data. Application of the low Reynolds number correc- 
tion yields good agreement with the data. 

Based on the comparison between predictions and data on a flat 
plate turbulent boundary layer, the following conclusions can be 
made : 

o The standard k- e model gives qualitatively good results. 
Significant improvements in mean velocity profiles are 
achieved by applying Chien’s low Reynolds number correc- 
tion to the k-e model (low Reynolds k-€) . 

o Algebraic stress model results are as good as the low 
Reynolds k-e model in regard to mean velocity profile, 

o The low Reynolds number correction is required for 
achieving good near-wall turbulent kinetic energy pro- 
files with both k-e and ASM models. 

o Individual fluctuating velocity components are reason- 
ably well correlated by ASM except in the viscous sub- 
layer, where significant improvements are obtained by 
applying the low Reynolds number correction. 
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ORIGINAL PAGE 18 
OF POOR QUALITY 


* a WATTS OATfl 
20 PARABOLIC 


A * WATTS DATA 



X=,846G 



A A 

-i — i — i — i — i — i— 1 

° f 

,mm r m 

— P- 

'~T Mm 


1 1 | 1 P J 

IQ 12 
U-VELOCITY 

14 16 IB 29 22 24 

(M/SEC) 

Q 

2 

4 

6 

8 10 12 
U-VELOCJTY 

14 t6 IB 20 22 24 

(H/SEC) 


A A WATTS DATA 



X«|.46S6 


1 “WATTS DATA 
20 PARABOLIC 


X=2.2743 














ORIGINAL PAGE 19} 
OF POOR QUALITY 



4 4 WATTS DATA 


*— 2D PARABOLIC 

X“.6408 


10 12 M IS IS 20 22 24 

U- VELOCITY (M/SEC) 


0 2 4 6 B 10 12 14 16 ' 8 20 22 24 

U-VELOCITY IH/SEC) 


4 4 WATTS DATA 

20 PARABOLIC 


4 4 WATTS DATA 
.04 J 2D PARABOLIC 


T , , r 

6 8 10 12 14 18 10 20 22 24 

U-VELOCI'v (M/SEC I 


-j , , r ..*M ftAj«-KAA-* T 1 r j. 

0 2 4 6 8 10 12 14 16 18 20 22 24 

U-VELOCl’Y (fi/SEC; 


4 4 WAT'S DATA 


2D PARASOL I C 

X= 1.4656 



i * 

•V j 

A A WATTS DATA 

A 

4 


A 

X* 1,8735 / 

i 


A 


a 


t> 

4*' , 

V ! 

A/ ! 

A/ 

A.' I 

... J 


A 2 4 B 8 10 12 14 16 18 20 22 24 0 2 4 6 8 10 12 14 16 18 26 22 24 

U- VELOCITY (M/SEC) U-VELOCI t t CM/SEC) 

Figure 6.1-3. Turbulent Boundary Layer Mean Velocity Profiles with 
the k-e Model and Chien’s Low Reynolds Number 
Correction (Low Reynolds Number k-e Model) . 










ORfQfNAL PAQg m 
OF POOR QUALITY 


* * MATTS BATA 


•— 3D PARABOLIC 

X-.A63* 


* AUATT5 DATA 
3D PARABOLIC 


1* 13 

1A IS 1* 33 33 3A 

• 

3 

4 

t 

8 1* 13 

1A 16 It 3« 

O-MELOCITY 

(R/SEC1 





U-UELOCITY 

(B/SEC) 


a A WATTS DATA 


'*• ED PARABOLIC 

X-.8AG6 



13 13 

JA 16 13 33 33 3A 

• 

3 

A 

G 

8 13 13 

1A 16 IS 33 33 

U-UELOCITY 

(B/SEC > 





U-UELOCITY 

(B/SEC) 


A A UATT5 DATA 


2D PARABOLIC 

X* 1.A65S 


A A WATTS data 


*— 3D PARABOLIC 

X-1.873S 


■ 2 A G 3 13 12 1A 15 II 31 S3 3A 

U-UELOCITY (A/SEC) 


SEAS 


13 13 1A 1G 18 38 33 3A 

U-UELOCITY (H/SECn 


Figure 6.1-4. Turbulent Boundary Layer Mean Velocity Profiles With 
the Algebraic Stress Model (ASM) . 
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6.2 Plane Couette Flow 


The plane Couette flow is a well defined flow for which anal- 
ytical solutions and good experimental measurements are available 
for evaluating the turbulence models. The test data selected for 

validating the models in this report were obtained by El Telbany 
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and Reynolds in a test setup shown schematically in Figure 
6.2-1. In this setup, the bottom wall was stationary and the top 
wall was moved at a velocity of 17.08 m/s, which corresponds to a 
Reynolds number of 12,640. The distance between the walls was 
44 mm. 

Computations for this flow were made with the standard k-e 
model and the algebraic stress model using a 2-D parabolic code. 
The standard k-€ model predictions and the data for mean velocity 
are shown in Figure 6.2-2. The agreement between data and predic- 
tions is very good. The non-dimensionalized Reynolds shear stress 
profile predicted by the standard k-e model is also in good agree- 
ment with the data, as seen in Figure 6.2-3. However, the profile 
of Reynolds stress normalized by the turbulence kinetic energy pre- 
dicted by standard k-e model is not in agreement with the data 
(Figure 6.2-4) * In the standard k-e model, in the regions where 
the shear stress is a constant, the turbulence kinetic energy is 
also a constant’. However, the data shows a gradual reduction in 
the k values away from the wall with the minimum value at the plane 
of symmetry. Consequently, the predicted uv/k profile is contant 
in the core of the flow, while the data shows a gradual increase in 
its value away from the wall. It is possible to match the predicted 
and measured values of uv/k at the plane of symmetry by increasing 
the turbulence model constant, C D , from 0.09 to 0.144. A signif- 
icant improvement is obtained in the uv/k profile, which is shown 
in Figure 6.2-5. However with C D = 0.144 used in the standard k-e 
model, the predicted mean velocity profile was not in agreement 
with the data, as seen in Figure 6.2-6. 
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The ASM predictions are illustrated in Figures 6.2-7 through 

6.2- 13. The comparison between ASM prediction and the data for 
mean axial velocity is shown in Figure 6.2~7. The agreement 
between the two is very good. The ASM prediction for uv, normal- 
ized by the wall shear stress is in good agreement with the data, as 
shown in Figure 6.2-8. However, when the Reynolds stress profiles 
are normalized by the local turbulence kinetic energy, shown in 
Figure 6.2-9, the predicted profile underestimates the values in 
the core of the flow. However, the ASM predictions for the center- 
line uv/k values are closer to the data than that predicted by 
standard k-e model. A very similar profile is obtained for the 
correlation coefficient, uv/(u'v'), which is shown in Figure 

6.2- 10. These two figures demonstrate that the ASM slightly over- 
estimates the turbulence kinetic energy components. 

The ASM prediction for the axial turbulence intensity, u', is 
shown in Figure 6.2-11. The predicted peak u* value is slightly 
smaller than the data. However, in the core of the flow, the ASM 
predictions are in good agreement with the data. The predicted and 
measured cross-stream turbulence intensity profiles are illus- 
trated in Figure 6.2-12. The predicted v* values are about 20 per- 
cent higher than the data. The predicted w* values are also higher 
than the data by about 15 percent in the region near the plane of 
symmetry as shown in Figure 6.2-13. 

The Couette flow calculations may be summarized as follows: 

o The standard k~e model predicts the mean velocity profile 
accurately, but underpredicts the centerline uv/k value 
by about 20 percent. 

o When the centerline uv/k value is matched with the data, 
(using C D = 0,144) , the predicted mean velocity profile 
is in poor agreement with data. 
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o 


The algebraic stress model correctly predicts the mean 
velocity profile and underpredicts the centerline uv/k, 
but the centerline uv/k values are in better agreement 
with the data than the standard k- e model. The basic 
reason for this deficiency is because of the overestima- 
tion of v* and w' by the ASM. Overall individual turbu- 
lence components are predicted well by the ASM. 
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Figure 6.2-1. Geometry of Plane Couette Flow. 
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Figure 6.2-3. Comparison Between k-e Model Predictions and 
Measured Shear Stress Profile (normalized by 
wall shear value). 
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6.3 Developing Flow in a Two-Dimensional Channel 


One of the simple flows considered for validating the k-« tur- 
bulence model was the developing flow in a two-dimensional channel. 
Analysis of the entrance flow problem provides a means of evalu- 
ating the accuracy of the numerical scheme. Detailed mean flow 

measurements in the entrance region of a parallel plate were made 
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by Emery and Gessner . The geometry of their test setup is shown 
in Figure 6.3-1. Predictions for this flow were obtained using a 
2-D elliptic code with the standard k-e model with 2200 grid nodes. 
Computations were performed until the total mass source error was 
less than 0.01 percent. Comparison between predicted and measured 
axial velocity variation along the centerline of the channel is 
shown in Figure 6.3-2. The difference between the two results is 
comparable to the measurement accuracy. Figure 6.3-3 illustrates 
the predicted and measured profiles of the axial velocity component 
at different axial stations* The agreement between data and pre- 
dictions is very good. 

The predicted and measured wall shear stress distributions are 
presented in Figure 6.3-4. The predictions and the measurements 
are within about 7 percent of the data, which is within the accur- 
acy of the wall sheai; stress measurements. The agreement between 
measured and predicted results demonstrates that the standard k-« 
model is sufficiently accurate for predicting mean flow field in a 
two-dimensional channel. 
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6.4 Developing Pipe Flow 

The next benchmark test case considered was the developing 

flow in a circular pipe. The test case selected for this problem 

68 

corresponds to the measurements by Barbin and Jones . The geome- 
try of their test setup is shown in Figure 6.4-1. 

The bulk flow velocity at the inlet in the test case was 
33.174 m/s and the Reynolds number/ based upon the bulk velocity, 
was 388,000. The mean velocity measurements were made using pilot 
tubes and the turbulence velocity components were measured using an 
x-wire probe. Computations for this case were made using a 2-D 
parabolic program. The computational domain extended from x = 0.3 
meters to x = 8.1 meters in the axial direction and from r = 0 to r 
= 0.1 meters in the radial direction. Along the axis of the tube, 
symmetry boundary conditions were specified, and along the pipe 
wall, standard wall functions were used to specify near-wall bound- 
ary conditions. In the computations, 100 grid nodal points were 
used in the radial direction. At x = 0.3 meters, the measured pro- 
files were used as initial profiles. Computations were made with 
standard k-« model and ASM. 

Comparison between standard k-€ model predictions and the data 
of Barbin and Jones for mean axial velocity is shown in Figure 
6.4-2. The mean velocity profiles are nondimensionalized by the 
average bulk velocity, = 33.17 m/s. The predicted mean velocity 
profiles are in very good agreement with the data. . 

The mean axial velocity profile comparison between the data 
and ASM predictions is presented in Figure 6.4-3. The ASM predic- 
tions are in good agreement with the measurements. The predicted 
and measured root mean square (RMS) axial velocity fluctuations, 
u', are illustrated in Figure 6.4-4. The ASM correctly predicts 
the axial turbulence intensity near the axis of the pipe. Near the 
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wall of the pipe, predicted peak values of u f are apparently higher 
than the data. However, the measurements near the wall do not have 
sufficient resolution. The boundary conditions near the wall were 
specified using wall functions, which have been shown (Paragraph 
6.1) to overestimate the peak turbulence kinetic energy in the case 
of a flat-plate boundary layer. By using an improved near-wall 
model, improvements in the peak turbulence intensity can be 
obtained. The comparison between the predicted and measured cir- 
cumferential turbulence intensity component, w', is illustrated in 
Figure 6.4-5. These profiles have characteristics very similar to 
the u' profiles. The w' peak values are slightly overestimated. 
The near-wall model deficiencies are responsible for the overesti- 
mation of the peak w* values. 


The results presented in this paragraph demonstrate that the 
k-e and ASM accurately predict the mean velocity profiles and that 
further improvements in turbulence structure and pressure drops can 
be achieved with an improved near-wall model. 
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Figure 6.4-1. Developing Pipe Flow Setup of Barbin and Jones 
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Figure 6.4-2 . Comparison Between k-e Model Predictions and Pro- 
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Figure 6.4-5. Comparison Between ASM Predictions and Measurements 

for Nondimensionalized Circumferential Turbulence 
Velocity Fluctuations. 



6.5 Fully Developed Pipe Flow 

The fully developed pipe flow is another case of simple flow 

where the turbulence models can be evaluated. The test measurement 
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selected for this case was that of Laufer , at a Reynolds number, 
R e , of 500,000. The geometry of the flow field is illustrated in 
Figure 6.5-1. Computations for this case were made with a 2-D 
parabolic program, starting with a plug flow profile at x = 0. The 
calculations were performed up to x = 10 meters, where fully devel- 
oped flow field was established. Predictions were obtained using 
the k-e model and the ASM with Chien's low Reynolds number correc- 
tion. A comparison between Laufer' s data and the k- e model predic- 
tions are shown in Figure 6.5-2. The agreement between data and 
predictions is very good. The predicted and measured turbulence 
kinetic energy profiles are shown in Figure 6.5-3. The standard 
k-e model predicts a higher value of peak turbulence kinetic energy 
near the wall compared to the data. Similarly, at the centerline, 
the k-e model predicts about 40 percent higher value for k than the 
data. 


The ASM prediction for mean velocity profile is shown in 
Figure 6.5-4. ' The predicted results and the data are in good 
agreement. The ASM prediction for turbulence kinetic energy is 
illustrated in Figure 6.5-5. The predicted peak as well as the 
centerline values of the turbulence kinetic energy are in good 
agreement with the data. The ASM predicts a faster decay of turbu- 
lence kinetic energy (k) , away from the wall, but the predicted 
variation of k in the core of the flow is slightly smaller than the 
measurements . 

The ASM predictions for axial turbulence intensity, a', and 
the data are presented in nondimens ional form in Figure 6.5-6. The 
predicted u' peak value near the wall is slightly smaller than the 
data. However, the agreement with data in the core of the flow is 
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very good. The predicted and measured v 1 profile is shown in 
Figure 6.5-7. The predicted profile is in good agreement with the 
data. However, the peak value of v' predicted by the model is 
slightly higher than the data. The predicted and measured w' pro- 
files, shown in Figure 6.5-8, are in good agreement in the entire 
flow field. At the axis of the pipe, the predicted w' value is 
slightly higher than the data. The comparison between data and 
predictions for the Reynolds shear stress, uv, is shown in Figure 
6.5-9. The data and predicted values are in excellent agreement. 

The low Reynolds number k-e model predicts the mean velocity 
profiles in a fully developed pipe flow accurately. It predicts a 
higher value of turbulence kinetic energy near the wall and at the 
centerline compared to the data. The ASM predicts the mean veloc- 
ity profile accurately and significantly improves the turbulence 
kinetic energy over the k— ^ model results • 
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Figure 6.5-3. Low Reynolds k- € Model Turbulence Kinetic 
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6.6 Two-Stream Mixing Laver 


Another benchmark test case selected from the data base for 
turbulence model validation is the flow in the mixing layer between 
two streams. Measurements for this case were made by Saiy and 
Peerless using a hot-wire probe and pitot tubes. A schematic of 
their flow test setup is shown in Figure 6.6-1. 

Computations for this case were made using the 2-D Parabolic 
Program. Since this flow field does not involve a wall boundary 
layer, low Reynolds number correction is not needed. Predictions 
for this case were obtained with the standard k- e model and the 
ASM. Initial conditions for these computations were applied at x = 
12.5 cm using measured data. A total of 100 cross-stream nodes 
were used in the computations. The nodes were closely distributed 
in the mixing region where gradients are higher and are sparsely 
spaced in the outer regions. For the test conditions, the veloc- 
ities of the two streams are; 

= 16.5 m/s; Uj = 38.37 m/s 

The predicted mean velocity and turbulence kinetic energy pro- 
files using the standard k- e model are presented in Figure 6.6-2. 

The predicted mean velocity profiles are in very good agreement 

\ 

with the measurements. The predicted turbulenbe kinetic energy 
(TKE) values are slightly smaller than the data. However, the 
width of the shear layer is correctly predicted. Overall agreement 
between k-e model predictions and data is good. 

Figure 6.6-3 shows the comparison between data and predictions 
obtained from the algebraic Reynolds stress model for the mean vel- 
ocity. The ASM predictions, similar to the k-e results, are in 
very good agreement with data. A typical comparison between data 
and predicted turbulence velocity components at x = 15 cm and x = 
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20 cm is presented in Figure 6.6-4. The u’ profiles are in very 
good agreement with the data. The predicted and measured v* data 
are in good agreement. A similar conclusion can be drawn for the w f 
component, with the exception that the predicted w' peak values are 
slightly larger than the data. 

The measured data of v' and w’ indicate that the peak values 
occur at y>o; i.e., they have shifted toward the low-velocity 
stream side of the mixing layer. The ASM model predicts them to lie 
along y = 0. 

Figure 6.6-5 illustrates the comparison between data and pre- 
dicted turbulence kinetic energy and shear stress (uv) profiles at 
x = 15 cm and x = 20 cm. The predicted turbulence k-e values are 
slightly smaller than the data. This is consistent with the 
results shown for the k-e model in Figure 6.6-2. The predicted uv 
profiles are in good agreement with datm except for a slight dis- 
crepancy at y = o. 

Major conclusions from the mixing layer work reported here 

are: 

o Both k- and ASM models give equally good mean velocity 
profiles as well as turbulent kinetic energy profiles. 
The measured peak values of the turbulent kinetic energy 
(I<E) are slightly higher than predictions; the ASM gives 
a little better correlation. 

o The ASM model gives good correlation for the fluctuating 
velocity components (u' f v', and w') as well as shear 
stress u ' v ’ ) . There is slight discrepancy in regard to 
the radial location of the v' and w' peaks. 
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Figure 6.6-2. Mixing Layer Mean Axial Velocity and Turbulent 
Kinetic Energy (TKE) Profiles Predicted by the 
Standard k-e Model. 
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6.7 Mixing of Coaxial Jets in Ambient Air 


Another benchmark test case selected from the data base for 
turbulence model evaluations is the flow in the near field uncon- 
fined mixing region of two coaxial jets. Measurements for the 

6 4 

selected test case were made by Champagne et al. using a hot-wire 

probe. A schematic of their test set-up is shown in Figure 6.7-1. 

For the test case studied, the ratio of outer to inner velocity at 

the nozzle exit was 5.0 with the area ratio, A^/A. = 2.94. 

o' 1 

Computations for this flow were performed using the 2-D para- 
bolic stream with the measured inlet velocity profile at x/D Q = 
1.16, where the maximum velocity, u was 18.29 m/s. The inlet 

lUclX 

kinetic energy profiles were obtained from measurements and a uni- 
form inlet length scale of 0.01 D q was prescribed. Computations 
were made with the standard k-e model and ASM. 

The predicted mean axial velocity profiles with the k-e model 
and data are presented in Figure 6.7-2 at x/D Q = 1.16, 2.14, 3.09, 
4.7, 6.05 and 8.02. The profiles shown at x/D Q = 1.16 are the ini- 
tial profiles used in the computation. Here YM2 represents the 
local half width of the jet. These results show that the data and 
k-e model predictions are in good agreement with each other. Fig- 
ure 6*7-3 show the comparison between data and ASM predictions for 
mean axial velocity. These profiles are in good agreement with the 
data, and a slight improvement over the k-e results can be seen. 

Figure 6.7-4 show the comparison between the data and ASM pre- 
dictions for u 1 . The predicted u' values are slightly higher than 
the data. However, the radial locations of the peak values are in 
good agreement with the data. Figures 6.7-5 presents the compari- 
son of the v' profiles. The predicted results and measurements are 
again in good agreement. The comparison between predicted uv and 
measured values are illustrated in Figure 6.7-6. These two results 
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are in very good agreement with each other up to x/D = 4.07. Beyond 
this station, the predicted uv values are slightly larger in magni- 
tude compared to the data. 

The k-€ and ASM predictions are in good agreement with mea- 
surements. Further improvements in ASM predictions can be achieved 
by fine tuning the empirical constants in the model. 
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Figure 6.7-5. Predicted and Measured Profiles of Fluctuating 
Radial Velocity Component (v 1 ). 





6.8 Free Circular Jet 


In partial support of the free swirling jet correlation and to 

further elucidate the coaxial jet mixing, a simpler case of free 

circular jet was selected for model validation. The benchmark case 
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selected for this flow was that of Wygnanski and Fiedler . They 
have reported accurate measurements of mean and turbulence velocity 
com^'onents using a hot-wire anemometer in a test setup shown sche- 
matically in Figure 6.8-1. Their jet diameter at the nozzle exit 
was 26.4 mm with a jet exit velocity of 51 m/s. Measurements were 
made at x/D = 40, 50, 60, 75, and 97.5. Computations for this case 
were made with a 2-D parabolic program using initial profile 
obtained from measurements at X/D = 40. Along the axis of the tube, 
symmetry conditions were applied. A total of 100 cross-stream 
points were used in the computations. 


The mean axial velocity profiles obtained from measurements 
and standard k-e model predictions are shown in Figure 6.8-2. The 
top left corner figure shows the initial profiles used in the com- 
putations. The predicted axial velocity results show a slower 

decay of centerline velocity than the measurements do. This may be 
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due to under-estimated diffusion rates. Launder has recommended 
modifying the turbulence model constants C D and C 2 for round jets 
in stagnant surroundings according to the relation 


C D = 0.09 - 0.04 f 


(141) 


where 


C 2 = 1.92- 0.067 f 



(142) 

(143) 


These modifications were used to predict the structure of 
Wygnanski and Fiedler's free jet. Comparison between the data and 
predictions for mean axial velocity are shown in Figure 6.8-3. The 
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predicted centerline velocity decay rate shown in this figure is 
smaller than the standard k-£ model results. 

The correction factor f is always positive and equations (141) and 
(142) would tend to reduce the magnitudes of C D and C 2 . Reduction 
of C D value would tend to reduce the eddy viscosity with attendent 
decrease in mixing rate. 

The present approach is to increase the C D value by the 
expression 


C D = 0.09 + 0.04 f (144) 

Furthermore, it was considered necessary to evaluate separating the 
effects of changing C c and C 2 « Application of equation (144) alone 
on the k-6 model will be denoted as the k-£l model and the use of 
equation (144) and (142) will be denoted as the k- £2 model in this 
report. 

The predicted axial velocity profiles using the k-ei model and 
the measurements are presented in Figure 6.8-4. The agreement 
between data and the k««€l model is excellent. Figure 6.8-5 allows 
the k-£2 model predictions for axial velocity. These results 
demonstrate that the k-€2 model tends to overestimate the mixing, 
which is responsible for the fast decay of the centerline velocity 
of the jet. 

The predicted mean velocity using the standard ASM is pre- 
sented in Figure 6.8-6. The ASM tends to slightly underestimate 
the mixing of the jet compared to the k-€l model, but is signifi- 
cantly better than the standard k-£ model. Comparison of the pre- 
dicted u* profiles and measurements, as shown in Figure 6.8-7, 
illustrates that u' is underpredicted up to x/D = 60 and beyond 
that station the agreement between data and ASM predictions is 
very good. This is a consequence of the underestimated mixing 
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rates in the model. Figure 6.8-8 shows the comparison of the pre- 
dicted and measured v* component at three axial stations. The 
agreement between these two is quite good. However, some disagree- 
ments can be seen close to the axis of the jet. The predicted w 1 
velocity profiles are in good agreement with the data, as seen in 
Figure 6.8-9. The predicted uv profiles are compared with the 
measurements in Figure 6.8-10. The uv values are initially under- 
predicted and are slightly overestimated at x/d = 75. 

For the case of the round free jet, the standard k~e model 
tends to underestimate the turbulent diffusion rates. Modifica- 
tions of the empirical constants are necessary to improve these 
results. The k-el model accurately predicts the mean velocity pro- 
files, while the k-€2 model tends to overestimate the jet center- 
line decay rate. The ASM shows a substantial improvement over the 
standard k-e model and no ad hoc modification of the empirical 
constants is necessary. The turbulence structure is well predicted 
by the ASM, and further refinement of the ASM is necessary to 
improve the quantitative predictions. 


Figure 6.8 
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Figure 6.8-3. Modified k-e Model Predictions with C D and C 2 
stants Given by Equations (141) and (142). 
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Figure 6.8-7. Comparison Between ASM Predictions and 
Measurements for Axial RMS Turbulence 
Velocity Fluctuations. 



RADIAL DISTANCE CM) RADIAL DISTANCE (M) 

0.100 0.200 0.300 0.400 0.500 0.000 0.080 0.160 0.240 0.320 0.400 


X/D - 60 X « 1.585 U 



0.00 0.50 1.00 1.50 2.00 

V-PRIME (M/S) 


ORIGINAL PAGE ESsj 
OF POOR QUALITY 


X/D - 75 


- 1.981 U 



V-PRIME (M/S) 


X/D - 97.5 X - 2,576 U 



0.25 0.50 0.75 1.00 

V-PRIME (M/S) 
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6.9 Flow Over a Heated Flat Plate 


The k~ e and the ASM seem to predict the momentum transport 
reasonably well for simple flows. It was deemed essential to eval- 
uate these models on the transport of scalar quantities such as 
temperature. The benchmark test case that was selected for this 
purpose was the flow over a heated flat plate. Measurements of 

mean and fluctuations of temperature were made for this case by 
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Charnay ' et al. A schematic of their setup is shown in Figure 

6.9-1. The flat plate was heated from the leading edge up to x = 

0.7 m and maintained at a uniform temperature of 313 °K. Beyond x = 

0.7 m, the wall temperature was abruptly changed and maintained at 

290 °K. The free stream t empe rature of air during the test was 

— . 2 

293°K. Measurements of T r T' and vT' were made at x = 0.7, 0.8, 
0„9, 1.05, and 1.4 m. 

Computations for this case were made using, the 2-D parabolic 
program with the initial conditions specified at x = 0.7 m from the 
measured Profiles. At x - 0.7 m, only the measured temperature pro- 
files were reported. The inlet velocity profiles were assumed to 

conform to the .law of the wall. The unknown wall mean stress was 

» 

calculated by assuming the temperature distribution to also follow 
a logarithmic law. The details of this calculation procedure are 
as follows: 

/ 

For a flat plate with constant free stream velocity and sur- 

u 

face temperature, the local Stanton number (- „ “ „ ) is given by 

•510 PooUooU . 


St 


- 0.2 

0.0287 Re 

X 

0.169 R - " 0 * 1 (13.2 Pr-10.16) +0.9 


where R a = p^U^ x/ju 
e x 


(145) 
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For the given free stream conditions, St was computed. The local 
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wall heat flux, q„ was obtained from 

% = Poo C p U (V T °o } St (146) 

The logarithmic law for temperature is given by Kays 2 " 1,2 as 

(147) 


(148) 

and y= £ 22 T kZ /no* 

From the prescribed temperature profile (T vs. y) , using equa- 
tion (148), the value of u* was computed. Knowing u*, the velocity 
profile was constructed from the law of the wall for mean velocity. 
This profile was used as the initial velocity profile. The turbu- 
lence kinetic energy was assumed to be a constant with k = 0.003 u . 
The length scale was assumed to be linear, 1 - k y up to y = 5 . 
Beyond that point, 1 was set equal to k 5. 

The boundary condition on the boundary layer edge was speci- 
fied through the computed entrainment rate. Along the wall 
boundaries, Chien’s low Reynolds number correction to the k-e model 
was applied. Across the boundary layer, a total of 100 nodes were 
distributed with the nodes closely spaced near the wall and further 
apart near the boundary layer edge. Computations were made using 

o k-e model with gradient transport model 

o ASM with gradient transport model 


t + = 2.195 In y + + 13.2 Pr - 5.66 

where, 

t + _ ( v t > u * 

(V^V 



o Algebraic scalar transport model (ASTM) 


The standard k-c model predictions with gradient scalar trans- 
port model are presented in Figures 6.9-2 through 6.9-4. In 
gradient transport assumption, the following expressions were used 
for the turbulent transports: 


where 



(150) 

(151) 

(152) 


Pr 


eff “ 


0.9 


Figure 6.9-2 shows the comparison between measured and pre- 
dicted mean temperature profiles across the boundary layer. In 
these figures, „the abscissa represents the temperature difference, 
(T~T Wq ) /(Too ~ T Wq ) , where T Wq is the wall temperature upstream of 
X = 0.7 m- (313 °K) . At X = 0.7, the nondimensional temperature 
profile would be similar to the velocity profile with a monotonic 
variation between zero at the wall to 1.0 at face stream edge. Just 
downstream of X = 0.7, the value of the nondimensional temperature 
at the wall jumps to 1.15. The mean temperature profiles gradually 
recover from a hot wall condition to a cold wall profile. The k-£ 
model predictions for temperature differences are smaller than the 
data as seen in Figure 6.9-2. In other words, the model under- 
estimates the heat transfer rate to the wall. 


The root mean square (RMS) value of the temperature fluctua- 
tions obtained from the k-e model are shown in Figure 6.9-3. In 
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these figures, the T* values are nondimensionalized by (T Wq -Too), 
(20°K in the present case) . The k-e model predicts high values of 
T' near the wall. At the outer edge of the boundary layer, the k-e 
model predictions are in reasonably good agreement with data. At 
x = 140 cm, the peak T' values tend to approach the measured 
values. This indicates that the gradient diffusion assumption is 
valid for equilibrium boundary layers. 

Figure 6.9-4 illustrates the comparison between data and k-e 
model predictions for the turbulent transport vT* . In these fig- 
ures, the quantity vT* is nondimensionalized by U (T Wq - T^) , with 
Tw q being the wall temperature upstream of x = 0.7 m. The k-e 
model underestimates the heat flux component v'T, especially in the 
region close to the wall. 

The predicted mean temperature profiles obtained from the ASM 
and gradient transport assumption are shown in Figure 6.9-5. These 
profiles are almost identical to those obtained from the k-e model, 
and the temperature differences are overestimated. The ASM pre- 
dictions for RMS temperature fluctuations are shown in Figure 
6.9-6. These profiles are also identical to these obtained from 
k-e model. A similar conclusion may be dravm for the vT' profiles 
obtained from ASM, as seen in Figure 6.9-7.. These figures illus- 
trate that the gradient transport model underestimated the heat 
flux, and the ASM does not signif icanlty improve the heat estima- 
tion. 


The predicted results using the ASTM are presented in Figure 
6.9-8 through 6.9-10. The ASTM uses the expressions given in Sec- 
tion 4.0 for the various turbulent transports. The ASTM predic- 
tions for mean temperature are shown in Figure 6.9-8. Comparison 
with the k-e model results (Figure 6.9-2) shows that the ASTM sig- 
nificantly improves the predictions for mean temperature, and at 
x = 140 cm, the predicted mean temperature profile agrees very well 
with the data. 
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The ASTM predictions for RMS temperature fluctuations are 

illustrated in Figure 6.9-9. The ASTM tends to overestimate the T* 

values near the boundary layer edge. In the near-wall region, some 

differences between the data and ASTM predictions are present. 

These differences are mainly due to the estimated turbulence struc- 
2 2 

ture, namely, u and v profiles. 

The ASTM predictions for the heat transport vT* are presented 
in Figure 6.9-10. These profiles are in good agreement with the 
data near the edge of the boundary layer. However, some dif- 
ferences exist in the near-wall region. These are due to the dif- 
ferences between estimated values and test conditions. The test 
results for the turbulence velocities were not reported. Further 
improvements in the ASTM predictions can be achieved if the turbu- 
lence structure predictions are refined. 

The results presented in this paragraph show that even for a 
simple flow case of boundary layer with sudden changes in wall 
temperature, the gradient transport assumption is not valid. The 
ASTM gives significantly improved predictions. Further improve- 
ments in the Reynolds stress predictions are needed to obtain 
quantitatively accurate results from ASTM. 


Comparisons of mean temperature (T-T Wq ) , the RMS temperature 
fluctuations ( T ' ) , and the heat transport (v^T") calculated using 
the various models ,can be made from the following figures: 

Figures 
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DELTA Y / DELTA 


ORKaSNAL PAGSjS 
. . - i t tTV 


o 



o 



o 



o 



VT/(UO(TWO-TINF)) x 10**4 


Figure 6.9-4. 
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Figure 6,9-6. ASM Predictions of RMS Temperature Profile. 
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Figure 6.9-9. ASTM Predictions of the RMS Temperature Profile 
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Figure 6.9-10. ASTM Predictions of (vT f ) Profiles. 
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6 .10 Plug Flow Reactor 


In reacting flows, the validation of the kinetic scheme is as 
important as the turbulence/chemistry interaction model. To vali- 
date the kinetic scheme, computations were made for the plug flow 
reactor shown schematically in Figure 6.10-1. 

Measurements in a plug flow reactor were conducted by Hautman, 
19 

et al., for lean, stoichiometric, and rich propane flames. These 
measurements were used to test the validity of the four-step scheme 
that has been proposed by Glassman and his associates based upon 
detailed species and temperature measurements under a well-con- 
trolled low-pressure and high inlet temperature environment. The 
Glassman four-step scheme has been incorporated into the Garrett 
Combustion Codes, both parabolic and elliptic. 

Computations were performed for lean, stoichiometric, and rich 
propane flames with both the two-step and the four-step schemes. 
Comparisons of these results with the measurements are shown in 
Figure 6.10-2 for the case of lean mixture. From Figure 6.10-1, it 
is clear that the four-step scheme is far superior to the two-step 
scheme in predicting the salient features of hydrocarbon combustion 
in the Princeton reactor. 

It should be noted that the four-step scheme as proposed by 
Glassman and his associates was based upon data from their plug 
flow reactor. This scheme probably represents a closer approxima- 
tion to actual hydrocarbon oxidation processes in a high tempera- 
ture environment than the simpler two-step scheme does. How four- 
step correlates other reacting flow situations, such as a laminar 
diffusion flame, premixed turbulent flames, and jet flames, is 
covered in Paragraphs 6.11, 6.12 and 6.13, respectively. 
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Figure 6.10-1. Princeton High Temperature Plug Flow Reactor. 
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Figure 6.10-2. Comparison of 2-Step and 4-Step Kinetic Scheme With 

Lean Propane Premixed Flame Data From High Tempera- 
ture Plug Flow Reactor . 









6.11 Laminar Diffusion Flame 


Another benchmark test case selected for validating the kin- 
etic schemes is the laminar diffusion flame. Measurements for a 
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laminar diffusion flame have been reported by Mitchell, et al. 
for the setup shown in Figure 6.11-1. This flow was computed with 
the 2-D elliptic CPM. Runs of this type are useful in the valida- 
tion of reaction mechanisms and establishing rate constant values 
since uncertainties due to turbulent interactions are absent. 

Comparisons between the measured and predicted species concen- 
tration, temperature, and velocity at different axial locations are 
shown in Figures 6.11-2 through 6.11-8, respectively, at three 
axial stations. The predictions were obtained with both kinetic 
schemes and include the influence of buoyancy and variable thermo- 
dynamic and transport properties. Overall, the agreement between 
the predictions and measurements is fairly good. 

Results with the two-step scheme for' the first axial station 
(x = 1.2 cm) are presented in Figures 6.11-2 through 6.11-4. The 
overall heat release rate as indicated by axial velocity (V) corre- 
lation is in gfrod agreement with data. Stable species profiles 
(e.g. C0 2 , H 2 0, 0 2 and unburned fuel) are also well correleted. 
The CO levels are predicted to be significantly lower^ than measure- 
ments by a factor of two to three. This is consistent with the two- 
step results on the plug flow reactor. Similar observation can be 
made for the comparison shown in Figures 6.11-5 and 6.11-6 at 
x = 2.4 cm. The temperature is slightly overpredicted (perhaps due 
to neglect of radiation losses in calculations) at x = 5 cm. The 
new CO prediction correlates well with the data. 

The overall agreement between the two-step predictions and 
data is reasonable. The slight discrepancies are due to the 
following: 
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(a) The presence of H 2 and intermediate hydrocarbons has been 
ignored in the two-step scheme. Tnis results in over- 
prediction of H 2 0, CH 4 , and temperature. Some of these 
differences can be overcome with an improved scheme. 

(b) Differential diffusion of species has been ignored. This 
will be significant for H 2 diffusion due to its low 
molecular weight. 

(c) Radiation was not considered in these computations. 

The four-step results with the rate constants suggested by 
Hautman et al. are presented in Figures 6.11-9 through 6.11-15. 
Except for the velocity profile at x = 1.2 cm, the scheme indicates 
good comparison with data as shown in Figure 6.11-9. The general 
shape of the temperature profile is well predicted as shown in Fig- 
ure 6.11-10. Near the flame centerline at x = 1.2 cm, the model 
underpredicts temperature levels by approximately 50 percent. The 
agreement improves further downstream, and- by x = 5 cm, the compar- 
ison is good. The fuel breakdown near the center is underpredicted 
by a factor of two as shown in Figure 6.11-11 for both x = 1.2 and 
2.4 cm. Similarly, for the initial portion of the flame, the C0 2 
levels near the center are considerably different from data. The 
same is true of the other stable species? for example, 0 2 and H 2 0 as 
shown in Figures 6.11-14 and 6.11-15. 

The four-step predictions in regard to CO are considerably 
better than the two-step results shown in Figures 6 . 11-4 through 
6.11-8. As shown in Figure 6.11-13, the CO peaks are similar in 
magnitude, as the data shows. Some improvement is desirable for 
the radial profile shape. Overall, the four-step correlates well 
with the data. Deficiencies are in regard to correlation in the 
initial portion of the flame centerline where the model under- 
predicts fuel breakdown and the levels of temperature, C0 2 , CO, 0 2 
and H 2 0. 
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Figure 6.11-1. Schematic of Laminar Diffusion Flame Setup Used By 

Mitchell. 
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Figure 6.11-3. 


Two-Step Predictions and Measurements for CH 
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Figure 6.11-4. Two-Step Predictions and Measurements for and 
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Figure 6.11-6. Predicted and Measured Species Profiles a 
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Figure 6.11-7. Comparison Between Predicted and Measured Axial 

Velocity and Temperature Profiles and Mitchell's 
Laminar Diffusion Flame at X = 5.0 cm. (Two- 
Step) 
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Figure 6.11-8. Predicted and Measured Species Profiles ai 
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Figure 6.11-9. 4-Step Scheme, Axial Velocity Profiles 
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Figure 6.11-10. 4-Step Scheme, Temperature Profiles 
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6 . 12 Turbulent Premixed Flame in a Rectangular Duct 


Another test case selected from the data Base for evaluating 

kinetic schemes is the reacting flow of premixed propane/air in a 

rectangular duct with a flame stabilizer as shown in Figure 6.12-1. 

137 138 

Measurements for this flow were made by Shipman, ' et al. 


Computations for this case were performed using the 2-D para- 
bolic program and standard k-e model, with the initial profiles 
obtained from measurements 0.0508 m downstream of the flame holder. 
The average inlet velocity was 18.288 m/s at a pressure of one 
atmosphere, having a turbulence intensity of 3 percent. Along the 

plane of symmetry, a zero radial gradient was specified for all the 

variables except V, which was set to zero. Along the outer radial 
boundary, wall function treatment was employed. Computations were 
performed with two-step and four-step kinetic schemes. On the two- 

X— « v*i W aw x. r-jj'-N .rll c?a j* n a •£" »• p +- o riAn e4*pnf^e ti»/a >* c> lie? o /~! f") n ra 

O OV^riCillC / lWU UJ.Xi.CJ. Ciiu dc uo vj. i U Uv nviv '-'i 

of them corresponds with the constants established in the Army Com- 
bustor Design Criteria Program, and the other set is for the PM/PV 


combustion. 
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Table 14 provides the values of the Arrhenius pre- 


expenents and the activation temperatures for the reaction steps, 


TABLE 14. TWO-STEP RATE CONSTANTS FOR GARRETT/AVLABS AND 
P REM I X ED/P RE VAPOR I ZING REACTION . 



REACTION STEP 1 

r 

REACTION STEP 2 

Pre- 

Exponent 

Activation 

Temperature 

Pre- 

Exponent 

Activation 

Temperature 

Design Criteria 
PM/PV 

3.3 x 10 14 
3 . 9xl0 9 

27.000 

18.000 

6 , 0x10® 
2 . 2xl0 8 

12,500 

12,500 


The predictions obtained using the Design Criteria constants 
are shown in Figures 6.12-2 through 6.12-8, and those using the 
PM/PV constants are presented in Figures 6,12-9 through 6.12-15, 


The predicted mean axial velocity using the Design Criteria 
constant and the measurements are shown in Figure 6.12-2 at five 
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axial stations downstream of the flame holder. The predicted axial 
velocity profiles show slower mixing rates compared to the measure- 
ments. The velocities near the plane of symmetry are considerably 
smaller in the predictions than in the data. 

Figure 6.12-3 illustrates the computed profiles of unburned 
fuel and data at five axial stations. The predicted unburned fuel 
mass fractions are higher than the measured values indicating 
slower fuel disappearance rate. Figure 6.12-4 illustrates the com- 
parison between predicted and measured CO mass fraction. The pre- 
dicted CO profiles are significantly lower than the data up to x = 
0.2032, partly due to the slow reaction rates. Beyond this sta- 
tion, the CO mass fractions are in reasonable agreement with data. 
However, the radial spread of CO profiles are underpredicted by the 
model. 

Since the reaction rates are underpredicted by the design cri- 
teria rate constants, the predicted temperatures (Figure 6.12-5) 
are also lower than the measurements. The radial spreading of tem- 
perature profile is also underpredicted by the model. 

The other derived variables such as 0 2 , C0 2 and H 2 0 are pre- 
sented in Figures 6.12-6 through 6.12-8, respectively. Due to the 
lower fuel consumption rate, predicted 0 2 profiles are higher than 
measurements. Similarly, the discrepancy between predictions and 
measured C0 2 can be explained. Apparent improvement in regard to 
H 2 0 profiles may be due to faster diffusion rates of this species 
compared to model assumptions of equal diffusivity for all species. 

In conclusion, the Army Combustor Design Criteria rate con- 
stants appear to underpredict fuel consumption rate and temperature 
profiles. On occasions the design criteria constants have seemed 
to overpredict reaction rates. More extensive validation is needed 
to establish the two-step rate constants for both diffusion and 
premixed flames. 
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The predicted results for mean velocity, obtained by using the 
PM/PV rate constants, are presented in Figure 6.12-9. A substan- 
tial improvement can be observed in the agreement between data and 
predictions, when compared to the results using Design Criteria 
rate constants (Figure 6.12-2). This is due to faster fuel con- 
sumption. Figure 6.12-10 illustrates the comparison of data and 
predictions for unburned fuel. These profiles are in much better 
agreement than the results obtained from the first set of rate con- 
stants. However, the radial diffusion rates are still 
underpredicted by the model. 

The comparison between predicted and measured C0 2 profiles are 
presented in Figure 6.12-11. Due to the improved convection rates, 
the reaction rates are expected to be higher, and hence the C0 2 
values are higher than the values obtained from the design criteria 
rate constants. The C0 2 values predicted from the PM/PV rate con- 
stants are still smaller than the measured values. 

The predicted and increased profiles for CO are illustrated in 
Figure 6.12-12. The predicted peak CO values are higher than the 
data and the predicted CO mass fraction profiles do not spread 
radially outwards as much as seen in the measurements. 

The predicted temperature distributions using the PM/PV rate 
constants and the measurements are presented in Figure 6.12-13. 
Since the reaction rates are faster, it is expected that the pre- 
dicted temperatures are also higher than those obtained using the 
design criteria rate constants. Overall these profiles are in good 
agreement with data. 

Figures 6.12-14 and 6.12-15 show comparisons between measured 
and predicted profiles of 0 2 and H 2 0, respectively. Although there 
is improvement over the design criteria constants in regard to 0 2 
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and C0 2 , the comparison is worse for H 2 0. This may be due to a 
number of reasons, such as: 

o Incorrect approximation of fuel breakdown by simple two- 
step 

o Neglect of H 2 as one of the intermediate products 

o Assumption of equal dif f usivities of all gaseous mole- 

cules 

o Turbulence/ chemistry interaction represented by a simple 
eddy breakup model 

Calculations were also performed using the four-step kinetic 
scheme outlined in Section 3.0. The rate constants used* in this 
computation were obtained from the report of Hautman, et al. 
These rate constants have given good comparison with plug flow 
reaction as shown in Figure 6.10-1. The Arrhenius pre-exponents 
and the activation temperatures for each of the four steps are 
given in Table 15. 


TABLE 15. RATE CONSTANTS FOR 4-STEP KINETIC SCHEME. 


REACTION 

STEP 

ARRHENIUS 
PRE-EXPONENT (K Q ) 

ACTIVATION 

TEMPERATURE 

EDDY BREAKUP 
CONSTANT (C R ) 

1 

2.0893 x 10 22 

24,800 

3.0 

2 

5.0117 x 10 19 

25,000 

3.0 

3 

3.9811 x 10 19 

20,000 

3.0 

4 

3.3113 x 10 18 

20,500 

3.0 



The predicted mean axial velocity profiles using the four-step 
kinetic scheme and the data are presented in Figure 6.12-16. The 
predicted velocities are significantly smaller than the measure- 
ments. Predicted unburned fuel profiles are shown in Figure 

6.12- 17 along with the data. Comparison with the two-step results 
(Figure 6.12-14) shows that there are no appreciable differences in 
the unburned fuel profiles between four-step and two-step. 

The four-step predictions for CO are presented in Figure 

6.12- 18. Due to the slow reaction rates in the four-step scheme, 
the predicted CO values are smaller compared to both the data and 
the two-step scheme (Figures 6.12-4 and 6.12-12). Consequently, as 
shown in Figure 6.12-19 the four-step predicted temperature pro- 
files are lower than the two-step (Figures 6.12-5 and 6.12-13) and 
the data. The other derived variables are similar including 0 0 , 

a 

CO 2 and shown in Figures 6.12-20 through 6.12-22. 

A number of reasons can be forwarded for delivering poor cor- 
relation with the four-step scheme. Numerical experimentation was 
made to demonstrate that the basic mechanism is valid and that 
future modifications to the approach will yield good comparison. 
Figures 6.12-23 through 6.12-29 present results with the first two 
reaction-step rate constants changed 

K 0l = 2.0893 x 10 24 Cr 2 = 6.0 

Ko 2 = 5.0117 x 1G 21 Cr 2 = 6.0 

Significant improvement in predictions can be seen, and one 
can therefore conclude that the basic four-step hydrocarbon oxida- 
tion mechanism is valid. 
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Figure 6.12-2. Predicted Velocity Profiles With Design Criteria 

Rate Constants. 


ORIGINAL PAGE IS 
OF POOR QUALITY 







Figure 6.12-3. Predicted Unburned Fuel Profiles With Design 

Criteria Rate Constants, 
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Figure 6.12-6. Predicted Mass Fraction With Design Criteria 

Rate Constants. 
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.12-7. Predicted C0 2 Profiles With Design Criteria Rate 
Constants. 
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Figure 6.12-9. Predicted Velocity Profiles With PM/PV Rate Con 

stants. 
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Figure 6.12-11. Predicted C0 0 Profiles With PM/PV Rate Constants 
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Figure 6.12-12. Predicted CO Profiles With PM/PV Rate Constants 
2 ' 
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Figure 6.12-13. Predicted Temperature Profiles With PM/PV Rate 

Constants. 
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.12-14. Predicted C >2 Profiles With PM/PV Rate Constants 
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Figure 6.12-18. Predicted CO Profiles with The 4-Step Scheme 
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Figure 

6.12-25. 

Modified 

4-Step CO Profiles 
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Figure 6.12-26. Modified 4-Step 
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Figure 6.12-29. 

Modified 

4-Step — H^O Profiles 





6.13 Free Methane Turbulent Jet Flame 


The free, turbulent, reacting methane jet investigated by 
133 

Hassan, et al., was modeled here using the PDF-partial eguilib- 

27 

rium approach of Bilger and Starner. Although several configura- 
tions were studied by Hassan, the detailed data was available only 
on the flame that has been investigated here. The test setup was a 
vertical free turbulent methane jet issuing from a 7.74 mm diameter 
pipe as shown schematically in Figure 6.13-1. The jet Reynolds 
number was 15,000, with a bulk velocity of 39.9 m/s. The flame was 
stabilized by a co-annular flow of 1 percent (by mass) of hydrogen. 
A discussion of the data and the results of the numerical model are 
presented below. 

From the initial (or tube exit) jet velocity, mass flow and 
assumed temperature of 59°F, mass and energy flux can be determined 
for comparison at downstream positions. The fuel was reported to 
be 94 percent methane. From the assumption of atmospheric pressure 
and the reported density, a molecular weight of 17.22 was deter- 
mined. If the remaining 6 percent was composed of nitrogen and 
propane (typical dry natural gas composition), then 52.4 percent of 
that fraction, being nitrogen would give the above molecular weight. 
From the reported mass flows of H 2 and fuel, the following mole 
fractions v/ere determined, 

Xpu e ] = 78.86%, ^CH4 = 74.13, ^N2 = 2.48, = 2.25, ^H2 = 21. 

The specific enthalpy at the jet exit is then ~4.154.10 6 and 

xg 

energy flux is -5858 — — . 
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At four axial positions (x/D = 36.2, 75, 113.7 and 204) radial 
profiles of temperature, mole fractions of CO, C0 2 and 0 2 were 
reported. Consider the reaction equation in mole fraction form, 

p C cx H cy N cn + A <°2 + sf N 2> ~ X C0 + X C0 2 + x 0 2 

+ x ch 4 + x ch 2 + x h 2 + X H 2 0 + x n 2 

Although dry measurements were made, the following analysis assumes 

that the reported results were based on the same reaction equation 

with the possible exception of CH 2 . Then, for the above reaction 

equation, there are four atom-balance equations and the sum of mole 

fractions identity for six unknowns (again exclude CH 2 ) ; therefore, 

an additional constraint is required. For this, the water/gas 
213 

equilibrium based on the measured local mean temperature and 
modified as described below, was chosen. Then at various radii in 
the cross section to. give an adequate definition, the reported data 
were interpolated using cubic splines and the above unknowns com- 
puted. Also the asymptotic end of each measured profile was deter- 
mined from cubic spline fitting (extrapolation) . If the above sys- 
tem of equations gave a negative mole fraction of methane (occur- 
ring in the high temperature and lean boundary regions) or CO was 
nonexistent, then the water/gas equilibrium constraint was dropped 
and the reaction equation solved (excluding both CH 2 and CH^) . 
Lastly, the outermost regions would not have any free hydrogen (as 
evidenced by computing negative mole fractions of same) . Then the 
carbon-balance constraint was dropped (taking the 0 2 measurements 
as being more accurate) , and the reaction equation solved for 
nitrogen, water, and stoichiometry. 

The water /gas equilibrium shift results of Mitchell, et 
129 

al., were applied to the equilibrium constant (determined from 
measured temperature) in the rich regions of the flame. In the 
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lean regions at less than 1400 °K, the reaction was assumed to be 
quenched {5) at 1400°K. 

With an assumed jet profile of [19], 



where the jet radius (r b ) is minimum of the C0 2 or 0 2 profile 
asymptote and the centerline velocity is determined by matching the 
energy flux. For this balance consider the following control vol- 
umes , 



Continuity gives ifi g - i“h^ + rh where ifi x is the entrainment, and 
energy gives 


fii 

h. + u i 
2 

r b 

" / 

2 

h + §- 

\ 

uprdr - 

1 

1 ^ 

f uprdr - ifii 



0 

_ _ 


0 


where , 

r fc - radius of thermal boundary as determined from 
the temperature profile 

^rad ” ra< 3iation heat transfer outside jet boundary 
to mass entrained 

p - local density 

itijL - initial jet mass flux 
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This procedure was attempted at all four measurement stations. 
For the first, x/D = 36,2, the energy flux matching gave u b = 26.85 
m/s. The other three stations have more enthalpy in the outer 
boundary region than matching will allow, indicating that temper- 
ature measurements there are high. 

The partial equilibrium computation%s were made with the full 
stoichiometry adjusted (i.e. amount of CH^ and N 2 ) to give the cor- 
rect enthalpy flux. The 2-D turbulent jet calculations were based 
on a jet diameter 1 mm larger than that reported (i.e., the OD of 
the H 2 stabilizer). Then, an assumed fully developed pipe flow 
velocity profile at the exit was set to give the correct jet mass 
flux. 214 The turbulence intensity profile was taken from developed 

n n c 

pipe flow results. The partial equilibrium computations were 

based on the specific thermodynamic of the JANAF tables using the 

2 16 

curve fits of Wakelyn and McLain, the three body recombination 

28 

kinetics of Jensen and Jones and the global hydrocarbon breakdown 
of Duterque, et al.** Additionally, the flame sheet approximation 
had a pyrolysis mixture fraction of 0.073 and 0.2 mass fraction of 
organic fuel converted to intermediate at the pyrolysis flame 
sheet. In the 2-D turbulent jet calculations, the initial turbu- 
lent length scales were 0.1 inner (jet) and 0.2 outer, both based 
on the jet diameter. Also the free stream or ambient air was given 
a velocity of f 0.25 m/s. 

The axial plots of centerline temperature, CO, C0 2 , and 0 2 are 
shown in Figure 6.13-2. The model-predicted centerline temperature 
profile agrees reasonably well in regard to temperature rise 
upstream and downstream of the flame tip at the center. The pre- 
dicted peak temperature level and its axial location are slightly dif- 
ferent from data. The initial CO buildup agrees well with data; 
but there is some discrepancy in the post-flame region. Similar 
conclusions can be made about the CO 2 profiles shown in Figure 
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6.13-2. For the initial portion of the flame (x < 0.75 m) , measure- 
ments indicate finite concentration of 0 2 at the centerline. Such 
a behavior cannot be predicted by diffusion name models including 
Bilger's model which, for a high fuel-rich region, sets 0 2 equal to 
zero . 


Figure 6.13-3 shows a comparison between measured and pre- 
dicted profiles of total fuel mass fraction, and unburned fuel pro- 
files are presented in Figure 6.13-4. In the initial portion of 
the flame, the conclusions are good, but farther down stream the 
fuel oxidation rate is faster than what data would indicate. This 
also results in higher centerline temperature levels as shown in 
Figure 6.13-5. From the model predictions of total fuel mass frac- 
tion at x/D = 113.7, it is concluded that the Bilger model is predict- 
ing a faster jet spreading rate. This causes faster decay of the 
centerline temperature in the post-flame region as shown in Figure 
6.13=5. 

Comparison between measured and predicted CO profiles (Figure 
6.13-6) show that, whereas the agreement is good up to x/D = 75, the 
post-flame region is not well correlated by the Bilger model. Sim- 
ilar conclusions can be made for the H 2 profiles as presented in 
Figure 6.13-7. For the 0 2 profiles, (Figure 6.13-8) up to x/D = 75, 
the model predictions are reasonable. But further downstream the 
model is predicting higher spreading rate than ' measurements sug- 
gest. Similar levels of correlations are obtained for CC> 2 profiles 
as shown in Figure 6.13-9. 
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Figure 6.13-2. Comparison' Between Bilger Model Predictions With 

Measured Centerline Profiles of Temperature, CO, 
CO ? and 0 2 for Hassan and Lockwood Methane Jet 
Flame . 
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Figure 6.13-3. Radial Profiles of Total Fuel (Bilger's Model) 
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Figure 6.13-4. Radial Profiles of Unburned Fuel Mass Fraction 
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Figure 6.13-8. . Radial Profiles of 0 2 (Bilger's Model) 
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